GROMACS 2026.0：NN势函数、GPU加速与AMBER/PLUMED完整支持

GROMACS 2026.0：NN势函数、GPU加速与AMBER/PLUMED完整支持摘要 GROMACS 2026.0于2025年1月19日发布，这是分子动力学模拟领域的一次重要突破。本文基于BioExcel Webinar #92的内容，为您梳理2026.0版本的核心亮点。六大核心更新：神经网络势函数接口：原生支持DeepMD、ANI等机器学习势模型，实现接近ab initio精度的经典MD速度 AMD GPU完整HIP后端：所有主要内核均支持AMD GPU，性能接近原生ROCm NVIDIA GPU自由能计算加速：FEP/TI的非键部分可在GPU上执行，性能提升10-30% AMBER力场完整验证：支持ff19SB、OL3等最新力场，与Amber软件完全兼容，用户可无缝迁移 PLUMED 2.9集成：增强采样功能更加稳定高效，长时间模拟不再崩溃 QM/MM稳定性改进：引入检查点机制，提高长时间模拟可靠性谁应该升级：GPU用户（AMD或NVIDIA）应立即升级以获得显著性能提升；需要高精度自由能计算或化学反应模拟的用户可以尝试NN势函数；AMBER用户现在可以无缝迁移到GROMACS，保留熟悉的力场参数；使用增强采样的用户将获得更稳定的PLUMED 2.9支持。视频信息来源：BioExcel Webinar #92 主讲人： Berk Hess（瑞典皇家理工学院 KTH） Lukas Müllender（瑞典皇家理工学院 KTH） Vedran Miletic（德国马普计算与数据设施）视频链接：https://www.bilibili.com/video/BV1Z3P4zeE4g，欢迎在bilibili关注『东山月光下』以观看视频，字幕已经上传！原始链接：What’s new in GROMACS 2026.0：https://www.youtube.com/watch?v=LUnOuUdTSwA 视频发布时间：2026年3月5日 GROMACS 2026.0发布时间：2025年1月19日核心亮点 1. 神经网络势函数接口这是2026版本最重磅的功能更新，它为GROMACS带来了机器学习势函数的原生支持，使得在经典分子动力学框架内运行接近ab initio精度的模拟成为可能。统一的接口设计：GROMACS 2026.0提供了通用的神经网络势函数接口，可以集成多种NN势模型，包括DeepMD、ANI、TorchANI等主流框架。这一接口的统一性意味着用户无需修改GROMACS源代码，只需提供训练好的模型文件即可使用。与GROMACS原生集成：接口直接使用GROMACS计算的pair list（邻接列表），避免了在NN模型内部重新计算非键相互作用，这是性能优化的关键。相比之下，许多外部NN势模型需要自己构建邻接关系，这在大型系统中会成为性能瓶颈。静电嵌入支持：接口支持QM/MM风格的静电嵌入方案，经典区域的电荷可以作为NN模型的输入，这使得NN模型可以感知周围经典原子的电场环境，从而实现更精确的QM/MM耦合模拟。这一特性对于研究化学反应、酶催化等需要量子力学精度的场景尤为重要。力反馈机制：NN模型计算的力可以作用于周围的经典原子，实现真正的双向耦合。这意味着NN区域和经典区域可以相互影响，而非简单的单向作用。对于蛋白质-配体复合物、溶剂化效应等研究，这一机制至关重要。工作流程：使用NN势函数的工作流程相对简单：首先需要准备训练好的NN模型文件（通常是PyTorch的.pt或.pth格式），然后在mdp文件中指定NN势函数模块并提供模型路径，GROMACS会自动加载模型并在运行时调用。 2. GPU性能飞跃 GROMACS 2026.0在GPU支持方面取得了革命性进展，不仅完善了对AMD GPU的支持，还在NVIDIA GPU上实现了自由能计算的加速。 AMD GPU完整HIP后端 2026.0提供了完整的HIP后端支持，使得GROMACS可以在AMD GPU上高效运行。HIP（HIP Interface for Portability）是AMD推出的GPU加速框架，旨在实现代码在AMD和NVIDIA GPU间的可移植性。完整的内核实现：相比之前的实验性版本，2026.0实现了所有主要内核的HIP后端，包括非键相互作用、PME长期静电、约束处理等。这意味着在AMD GPU上运行GROMACS不再需要功能妥协，可以获得与NVIDIA GPU相当的完整功能体验。性能接近原生ROCm：根据官方测试，HIP后端的性能接近AMD原生ROCm优化代码，在某些场景下甚至可以达到90%以上的性能。这一性能水平已经足以满足大多数生产环境的需求。严格的测试验证：HIP后端经过了系统的单元测试和集成测试，不仅由GROMACS团队在标准测试基础设施上验证，还由AMD开发人员进行了独立测试。目前HIP后端的性能已达到相当成熟的水平，可以放心用于生产环境。 NVIDIA GPU自由能计算加速 GROMACS 2026.0将自由能计算内核移植到了CUDA GPU上，这是继PME和键长约束之后的又一个重要GPU加速模块。非键自由能内核GPU实现：自由能微扰（FEP）和热力学积分（TI）等方法的非键相互作用部分现在可以在GPU上执行。这包括Lennard-Jones势、库仑相互作用等的自由能微扰项。之前这些计算必须在CPU上完成，成为性能瓶颈。 CPU-GPU异步执行：GPU和CPU可以并行工作，GPU计算非键自由能贡献的同时，CPU可以处理其他任务。这种异步执行模式在GPU很快、CPU相对较慢的配置下性能提升尤为显著。适用场景：自由能GPU加速在以下场景下效果最佳：当你有快速的GPU和相对较慢的CPU，或者你扰动了系统的很大一部分原子（如大分子配体的结合）。在典型的小分子自由能计算中，性能提升可达10-30%。为什么之前没做：很多人可能会问，为什么GROMACS没有早点实现这个功能？原因是在很多情况下，CPU在GPU计算时是空闲的，将自由能计算放到GPU上并不能提升总体性能。但随着GPU速度越来越快，CPU-GPU性能差距扩大，GPU加速自由能计算变得有意义了。多GPU性能优化对于拥有多GPU的高端系统，2026.0引入了GPU-direct通信和多rank PME等重要优化。 GPU-direct通信：在多GPU模拟中，GPU之间的数据传输（如PME网格交换）现在可以通过GPU-direct技术直接进行，无需经过CPU内存。这大大降低了通信延迟，提高了带宽利用率。多rank PME在GPU上并行：PME（Particle Mesh Ewald）长期静电计算的多个rank可以在GPU上并行执行，充分利用多GPU的计算资源。性能提升：在标准测试中，多GPU优化带来了5%的性能提升。虽然数字看起来不大，但在长时间模拟中累积下来仍然是显著的提升，特别是对于大规模生产模拟而言。 3. AMBER力场完整集成与验证 GROMACS 2026.0对AMBER力场的支持进行了系统性的改进和验证，确保与Amber最新版本的兼容性。包含最新AMBER力场：2026.0支持ff19SB蛋白质力场、OL3 RNA力场等AMBER最新版力场。这些力场代表了AMBER力场家族的最新进展，在蛋白质和RNA的模拟精度上有显著提升。完整的验证流程：GROMACS团队对新版AMBER力场进行了系统的测试和验证，包括小分子、蛋白质、核酸等多种测试体系。验证工作不仅由GROMACS团队完成，还得到了AMBER开发团队的确认，确保与Amber软件的计算结果一致。参数兼容性保证：用户现在可以放心地将在Amber中构建的模型迁移到GROMACS，不用担心力场参数的差异。这对于需要同时使用两个软件的用户（例如在Amber中做参数化，在GROMACS中做生产模拟）来说是一个重大利好。 4. PLUMED增强采样集成更新 PLUMED是分子动力学增强采样的核心插件之一，GROMACS 2026.0更新了对最新PLUMED版本的支持。更新至PLUMED 2.9：集成了PLUMED 2.9版本，这是PLUMED项目的最新稳定版本。PLUMED 2.9带来了许多新功能和性能优化，包括新的偏置势方法、改进的元动力学算法等。不是2.10.0吗？改进的集成接口：GROMACS与PLUMED之间的接口更加稳定和高效，降低了崩溃和内存泄漏的风险。这对于长时间增强采样模拟尤为重要，因为这类模拟通常需要运行数天甚至数周。支持更多模块：更新后的接口支持更多PLUMED模块和势函数，包括用于研究蛋白质折叠、配体结合、相变等过程的专用模块。用户可以更灵活地设计增强采样策略。 5. 运行时性能监控指标 GROMACS 2026.0在日志文件末尾添加了新的性能指标，帮助用户更好地评估和优化模拟性能。每步毫秒数（ms/step）：显示每一步MD模拟所需的毫秒数，这是最直观的性能指标。通过监控ms/step，用户可以快速判断模拟是否达到预期性能，以及是否存在性能瓶颈。每秒百万原子步数（$10^6$ atoms × steps/s）：这是一个归一化的性能指标，综合考虑了体系大小和模拟速度，便于在不同大小的系统之间比较性能。数值越高说明模拟效率越高。这些指标在日志文件末尾自动输出，用户无需手动计算，大大简化了性能评估工作。特别是在尝试不同参数组合时，这些指标可以帮助快速找到最优配置。 6. QM/MM稳定性改进对于使用QM/MM方法的用户，GROMACS 2026.0引入了一个看似微小但影响重大的改进：QM中心定位的检查点（checkpointing）功能。问题背景：在之前的版本中，如果QM中心在模拟过程中偏离初始位置太远，系统可能会变得不稳定，甚至导致模拟崩溃。这是因为QM区域的定位信息没有被保存和恢复。检查点机制：2026.0实现了QM中心定位的检查点功能，当写入检查点文件时，QM中心的坐标和定位信息会被保存。从检查点恢复模拟时，这些信息会被正确恢复，确保模拟的连续性和稳定性。实际影响：对于长时间QM/MM模拟或需要频繁重启模拟的用户，这一改进大大提高了模拟的可靠性。你不再需要担心因为检查点问题导致模拟失败，这在生产环境中是一个重要的稳定性保证。版本号规则解读从2026版本开始，GROMACS采用全新的版本号规则，这一变化旨在让版本号更加直观和一致。主版本号：年份（如2026）表示主要功能发布版本。每年通常会发布一个主版本，包含新功能、性能优化等重要更新。次版本号：bug修复版本（如2026.1、2026.2）只包含错误修复和文档改进，不添加任何新功能。这确保了次版本升级的稳定性，用户可以放心升级而不用担心功能变化带来的兼容性问题。升级建议：建议始终使用最新的次版本号，因为bug修复可能解决你遇到的问题，而且不会破坏现有工作流程。例如，如果你使用2026.0，遇到bug后应该升级到2026.1或更高版本，而不是停留在旧版本。适用场景与实用建议神经网络势函数适合这些场景需要ab initio精度但经典MD速度的研究：例如研究化学反应机理、酶催化过程、电子结构敏感的性质等。NN势函数可以提供接近DFT精度的能量和力，但计算成本接近经典力场。复杂化学反应研究：NN势函数可以处理键断裂和形成过程，这是传统经典力场无法做到的。例如研究蛋白质折叠过程中的二硫键形成、小分子在酶活性中心的反应等。高精度自由能计算：使用NN势函数计算结合自由能、溶剂化自由能等，可以获得更可靠的结果。对于药物设计领域的用户，这意味着更准确的亲和力预测。 QM/MM耦合模拟：NN势函数可以替代传统的QM区域，提供更低成本但保持足够精度的量子力学描述。特别适合大型生物分子的QM/MM模拟。 GPU加速适合这些场景大规模体系（>10万原子）：例如膜蛋白-脂质双分子层体系、核糖体等大分子复合物、病毒衣壳等。GPU加速可以大幅提升这些体系的模拟速度。长时间尺度模拟（微秒级）：GPU加速使得微秒级模拟在合理时间内完成成为可能。例如研究蛋白质构象变化、膜蛋白-配体结合动力学等需要长时间采样的过程。多GPU并行计算：对于拥有多GPU的工作站或集群，2026.0的多GPU优化可以充分利用硬件资源，获得接近线性的性能提升。自由能计算：自由能微扰、热力学积分等计算密集型方法在GPU上的加速尤其明显。对于需要计算多个配体的结合自由能的药物设计项目，GPU加速可以节省大量计算时间。参考资源 GROMACS官网：https://www.gromacs.org/ BioExcel网站：https://bioexcel.eu/ 视频链接：https://www.youtube.com/watch?v=LUnOuUdTSwA GROMACS手册：https://manual.gromacs.org/ 论坛讨论：https://gromacs.bioexcel.eu/ 字幕翻译与整理：东山月光下（B站）。本文基于BioExcel Webinar #92的字幕整理而成

Molecular Dynamics · 2026-03-06

Amber ff19SB高温MD模拟的水模型选择、系综设置与金属离子参数

Amber ff19SB高温MD模拟的水模型选择、系综设置与金属离子参数搜到的资料不多，结合了AI整理和推断，如有错误恳请指出[合十][合十]。摘要在高温分子动力学模拟和金属离子体系建模中，水模型选择、系综设置和离子参数配套共同决定模拟结果的可靠性。本文系统性地梳理了 OPC 与 OPC3 的适用边界、450 K 高温构象采样的系综选择逻辑，以及高价金属离子的 12-6-4 模型参数化与验证。对于水模型选择，ff19SB 论文在已测试水模型中推荐与 OPC 组合（未评测 OPC3）；独立基准研究显示 OPC 在宽温区密度–温度曲线和热膨胀系数上整体优于 OPC3。对于 450 K 构象探索，推荐使用 300 K NPT 确定密度后进行 NVT 高温采样，最终回到 300 K NPT 重新平衡[3]。对于三价/四价金属离子，传统 12-6 模型无法同时重现水化自由能（HFE）与离子–氧距离（IOD），误差可达 ±100 kcal/mol（HFE）和 ±0.1 Å（IOD），必须使用包含 $C_4$ 项的 12-6-4 模型（误差分别在 2 kcal/mol 与 0.01 Å 以内）。在超氧化物还原酶（$\ce{Fe^{3+}}$ + OPC）的验证中，图8 和图9 共同证明：12-6-4 模型在保留配位球结构方面显著优于 12-6 模型，且优化 IOD 的 12-6 参数集在配位几何稳定性上也优于 12-6 HFE 参数集[5]。更换水模型时必须同步配套对应的离子参数，否则可能导致系统性偏差。核心结论水模型优先级：ff19SB 原论文在已测试的显式水模型中推荐 ff19SB + OPC，且未评测 OPC3；若受限必须使用三点水，可选择 OPC3 作为折中方案[4] 高温性能判断：基准研究显示 OPC 在宽温区密度–温度曲线和热膨胀系数上整体优于 OPC3；12-6 模型下 OPC3 的 IOD–HFE 曲线最接近实验目标点，但仍有系统性误差[1][2][5] 构象采样策略：450 K 用于初始构象探索时，建议以 300 K NPT 的体积进入 NVT 高温采样，最终结论以 300 K NPT 的再平衡与生产采样为准[3] 离子参数配套：更换水模型后必须同步更新对应的离子 Lennard-Jones 参数；对于三价/四价金属离子，优先采用 12-6-4 模型，其定量优势在图5部分详细说明[5] 12-6-4 在蛋白体系中的验证：在超氧化物还原酶（$\ce{Fe^{3+}}$ + OPC）的验证中，图8 和图9 共同证明12-6-4在保留配位球结构方面显著优于12-6；且优化IOD比优化HFE更重要，12-6 IOD参数集的配位几何稳定性远优于12-6 HFE参数集[5] 物理机制：OPC 的 M-site 有助于更好拟合高阶多极矩，从而改善氢键网络与温度依赖性质[1][2] 背景高温分子动力学模拟（如 450 K 退火或加速采样）在蛋白质构象探索和增强采样中广泛应用。然而，高温条件下的水模型选择往往被研究者忽视，导致模拟结果可能引入不必要的系统偏差。水模型作为 MD 模拟中占比最大的组分（通常占体系原子数的 80% 以上），其性质对体系的动力学行为、热力学响应和溶剂化结构具有决定性影响。在常温（300 K）下，大多数主流水模型（TIP3P、OPC、OPC3 等）都能给出合理的结果。但在高温或宽温区研究中，不同水模型对温度依赖性质（如密度随温度的变化、热膨胀系数、介电常数等）的拟合能力差异显著。当前存在一个关键的知识缺口：当研究者需要使用 Amber ff19SB 这一代高精度蛋白力场进行高温 MD 模拟时，应该选择 OPC 还是 OPC3 水模型？两者在 450 K 下的性能有何差异？在 NVT 和 NPT 系综之间应该如何选择？这些选择背后的物理机制是什么？水模型选择 ff19SB 水模型选择：OPC 还是 OPC3？在设计高温 MD 模拟方案时，第一个需要明确的问题是：ff19SB 力场应该搭配哪个水模型？ ff19SB 的水模型兼容性 ff19SB 力场以氨基酸特异的 CMAP 修正主链 $\phi/\psi$ 能量面，共拟合 16 组 CMAP（$24 \times 24$ 网格），训练目标为溶液相 QM 能量面，因此不依赖于某一个固定水模型。从兼容性角度，ff19SB 可以与 OPC、OPC3、TIP3P 等多种水模型组合使用。 ff19SB 原论文仅比较了 OPC 与 TIP3P 并推荐在已测试的显式水模型中使用 OPC，同时强调 ff19SB 并未用 OPC 拟合，水模型仍可能是限制因素，未来其他水模型不排除更好[4]。需要说明的是，OPC3 并未包含在 ff19SB 原论文的评测范围内，本文关于 OPC3 的讨论主要来自水模型基准研究。 http://archive.ambermd.org/202303/0144.html 里提到[6] Hi Vlad, Yes we have done some tests using opc3, nothing published yet. For peptides the match to experiment degrades a little compared to opc, but better than tip3p. I don’t have more specifics since I am at the ACS meeting this week. Carlos OPC vs OPC3：本质区别 OPC（Optimal Point Charge water）与 OPC3（Optimal Point Charge 3-point water）是同一研究团队开发的两种水模型，它们的本质区别在于点位（sites）布置和电荷分布方式：特性 OPC OPC3 点位类型 4-point 模型 3-point 模型电荷布置除了两个 H 和 O 以外，还有一个无质量的负电荷点（M-site）偏离氧原子中心，O上无电荷所有电荷都放在 O/H 原子上电荷参数 q=0.6791 e[2] q=0.447585 e[1] 几何参数 l=0.8724 Å，$z_1$=0.1594 Å，θ=103.6°[2] l=0.97888 Å，θ=109.47°[1] LJ 参数 $\sigma_\mathrm{LJ}$=3.16655 Å，$\varepsilon_\mathrm{LJ}$=0.89036 kJ/mol[2] $\sigma_\mathrm{LJ}$=3.17427 Å，$\varepsilon_\mathrm{LJ}$=0.68369 kJ/mol[1] 设计理念类似 TIP4P 的思路，通过 M-site 更准确地拟合水分子的静电分布与氢键网络在 3 点刚性水模型的精度上限约束下做的最优拟合拟合目标优化整体水性质和溶质–水相互作用在 3 点模型框架下达到最佳拟合注：$z_1$ 表示负电荷虚拟点（M-site）相对氧原子沿水分子对称轴的位移，OPC3 为三点模型因此不适用。[1][2] 两者的共同点是以电荷分布为核心进行优化。OPC 的构建采用对 $\mu$–$Q_T$ 空间的系统搜索，仅保留对称性约束，以优化液相电静特征；OPC3 在相同思路下将模型压缩为三点形式，以获得更高的计算效率[1][2] 从物理意义上理解，OPC 的 M-site 相当于在氧原子附近增加了一个额外的“虚拟电荷点”，使得模型能够更准确地再现水分子的高阶多极矩（quadrupole moment），从而改善对氢键网络和溶剂化结构的描述。这里的 $\mu$ 表示水分子偶极矩，$Q_T$ 表示四极矩的迹。OPC 论文定义了一个质量评分，用多项体相性质与水化自由能的综合误差来衡量模型在 $\mu$–$Q_T$ 空间的优劣，得分越高表示越接近目标性质[2]。图1：OPC 的 $\mu$–$Q_T$ 质量评分图（原文 Figure 3）[2] 该图展示了在 $\mu$–$Q_T$ 空间中的模型质量分布，OPC 位于高质量区域，说明其电静多极矩选择更接近液相最优区间[2]。精度 vs 速度/兼容性 OPC 和 OPC3 的选择本质上是在模拟精度与计算通用性之间做权衡： OPC 的优势：在整体水性质、溶质–水静电相互作用、氢键网络的再现上通常更准确。但 4 点模型在某些 MD 引擎或工作流中会稍麻烦或略慢（如 GPU 加速路径对 4 点水的优化程度可能不如 3 点水）。 OPC3 的优势：通常更快、更“通用”（3 点水对很多程序/加速路径更友好），但就水本身的综合性质拟合而言一般不如 OPC。社区实践经验基于原论文结论与常见实践，若不受 3 点水限制，优先使用 OPC；若必须使用 3 点水，再以 OPC3 作为替代。 ff19SB + OPC 的实验验证：图11：CLN025 蛋白的主链 RMSD 随时间变化（Maier et al., JCTC 2020, Figure 11）[4] 该图展示了在 CLN025（一种快速折叠的 β-hairpin 蛋白）的模拟中，三种力场+水模型组合的性能：从天然结构（nat）与完全伸展结构（ext）出发，各 4 条轨迹，共 8 次独立模拟；300 K 进行，总时长约 172 μs 性能对比： ff19SB + OPC（蓝色）：能够可逆地折叠到天然结构，native population = 50 ± 17% ff14SB + TIP3P（红色）：native population = 75 ± 23% ff14SB + OPC（黄色）：native population = 33 ± 19% 关键发现：折叠可逆性：4 次 nat 与 4 次 ext 轨迹均回到天然结构，说明该组合稳定可靠组合匹配性：ff14SB + OPC 的 native population 低于 ff14SB + TIP3P，提示 OPC 与 ff14SB 的协同不足协同优势：ff19SB 并未专门拟合 OPC，但与 TIP3P 对比时 OPC 在折叠动力学与构象平衡上更好[4] 这个实验数据支持 ff19SB + OPC 作为推荐组合的结论，特别是在蛋白折叠、构象平衡等应用中[4]。一个实用的经验法则：默认（蛋白折叠/构象平衡/IDP 等）：ff19SB + OPC 必须 3 点水（例如某些代码路径、极限性能、或你工作流只能稳定支持 3 点）：用 OPC3，并确保离子参数选择合理/一致高温下的性能差异：OPC 还是 OPC3 更好？高温（450 K）是水模型性能差异被放大的场景。当温度升高，水分子的动能增加、氢键网络减弱、密度下降，不同水模型对温度依赖性质的拟合能力差异会显著影响模拟结果的可靠性。纯水基准测试：宽温区对比多项研究已经系统对比了 OPC 和 OPC3 在宽温区（270–650 K）的表现： OPC3 相关论文（Izadi & Onufriev, 2016）：直接对比了 OPC vs OPC3 的密度–温度曲线，作者明确指出：[1] 4-point OPC 在宽温区密度的温度依赖上比 3-point OPC3 更准确给出了一个关键的派生量：OPC3 的热膨胀系数偏差（约 $67.9\%$）远大于 OPC（约 $5\%$）文中指出 OPC3 在三点模型中显著优于 TIP3P/SPC/E，并认为实用三点刚性非极化模型已接近精度上限 2024 年三点水模型的大规模对比（11 个刚性三点水模型）系统评估了液–汽共存、临界点与自发气化等高温行为：[3] 给出各模型的 $T_\mathrm{C}$、$T_\mathrm{MD}$ 与 $T_\mathrm{evap}$，$T_\mathrm{evap}$ 范围约为 $520$–$620~\mathrm{K}$，并明确指出 $T_\mathrm{evap}$ 不是沸点该研究仅覆盖三点模型（包含 OPC3），不包含四点 OPC，因此不能据此得出 “OPC3 优于 OPC” 的结论 OPC 原始论文强调：OPC 通过优化点电荷分布来逼近液相电静特征，体相性质平均相对误差约 $0.76\%$，并且在宽温区保持与实验接近；同时小分子水化自由能的 RMS 误差可做到 $<1~\mathrm{kcal/mol}$[2]。高温性能差异从何而来？ OPC vs OPC3 在高温下的性能差异，核心来自电荷点位布置的不同： OPC（4-point，带 M-site）：负电荷不锁死在氧原子上，而是分布在 M-site → 能更好复现高阶多极矩，从而改善氢键网络与温度依赖性质 OPC3（3-point）：负电荷必须在氧上 → 多极矩表达受限，作者明确指出这会拖累密度温度依赖与热膨胀等指标[1] OPC3 论文给出了两者的多极矩差异：OPC 的 $\mu = 2.48~\mathrm{D}$、$Q_T = 2.3~\mathrm{D\cdot Å}$，而 OPC3 的 $\mu = 2.43~\mathrm{D}$、$Q_T = 2.06~\mathrm{D\cdot Å}$[1][2]。 OPC 的负电荷可偏离氧原子以更好兼顾高阶多极矩；OPC3 负电荷固定在氧上，导致高阶多极矩拟合受限。直接回答“高温下谁更好？” 如果你说的“高温”是指温度高于 350 K 甚至更高并且你关心温度依赖的体相水性质：倾向选择 OPC 如果你受限于 3 点水（性能/引擎/工作流），OPC3 是可接受的折中方案，但要接受它在密度–温度曲线/热膨胀上偏差更大。 450 K 构象采样：NVT 还是 NPT？当你的研究目标是 450 K 下进行蛋白质构象采样（如高温退火、加速跨越能垒），系综的选择（NVT vs NPT）和体积/密度的设定策略会直接影响采样效率和结果可靠性。 NVT vs NPT：物理意义的本质区别首先需要明确 NVT 和 NPT 系综在高温下的物理含义： NVT（等温等容）：固定体积，温度耦和到热浴。体系密度被锁死，不会因温度升高而膨胀。 NPT（等温等压）：固定压力（通常 $1~\mathrm{bar}$），体积可以自由调整。体系会根据温度自动调整到平衡密度。在 $450~\mathrm{K}$、$1~\mathrm{bar}$ 的条件下，液态水处于超热液体区域。对 11 种刚性三点水模型的系统研究表明，NPT 下存在模型相关的自发气化温度 $T_\mathrm{evap}$，且 $T_\mathrm{evap}$ 并不等于沸点。该研究给出的 $T_\mathrm{evap}$ 范围约为 $520$–$620~\mathrm{K}$，其中 $T_\mathrm{evap}$ of OPC3 为 $593.7 \pm 1.2~\mathrm{K}$（C-rescale barostat）[3]。因此，450 K 低于 $T_\mathrm{evap}$，体系在 NPT 下仍可能保持液相，但密度会明显下降，并对 barostat 与升温速率更敏感。若继续升温接近 $T_\mathrm{evap}$，则可能出现空泡、密度骤降、体积迅速增大的“自发气化”现象。你关心的问题类型选择 NVT 还是 NPT，取决于你的研究目标： 1) 只是要一个稳定溶剂环境（重点关注蛋白高温退火/加速采样） ✅ NVT 是合理选择。OPC3 可以用（或 OPC，如果你能用 4-point）。作为三点模型，OPC3 在温度依赖的体相性质上精度有限，但用于“稳定溶剂环境”的需求通常足够。在这种用途里，决定能否稳定运行的往往不是水模型，而是：初始密度是否合理（NVT 下密度不会自动纠正）约束/时间步/恒温器设置是否稳定一个常见参照是温度‑REMD：多数 REMD 实现会在 NVT 下运行多个 replica，在 Amber 这类力场工作流中也很常见；Amber 早期 REMD 只支持 NVT，后续才扩展到 NPT‑REMD[7][8]。因此，把高温 NVT 当作构象探索的工具是合理的，但最终统计仍应回到常温 NPT 的再平衡与生产采样。如果你只需要“稳定液相环境”，核心问题是 $450~\mathrm{K}$ 是否低于 $T_\mathrm{evap}$。三点水模型的大规模对比研究给出 OPC3 的 $T_\mathrm{evap}=593.7 \pm 1.2~\mathrm{K}$，明显高于 $450~\mathrm{K}$，因此在 $450~\mathrm{K}$ NVT 下使用 OPC3 作为稳定溶剂环境是合理的[3]。需要强调的是，高温轨迹只用于初始构象探索，最终统计应回到 $300~\mathrm{K}$ NPT 重新平衡与生产采样。若进行高温 NPT 预平衡，建议采用 C-rescale 并先在中间温度预平衡密度。 2) 你要在 450 K 下比较水的热力学/界面性质（密度-温度曲线、热膨胀、表面张力等） ⚠️ 需要谨慎：OPC3 论文认为实用三点刚性非极化模型已接近精度上限；相比之下 OPC（4-point）在密度温度依赖与热膨胀上通常更贴近实验[1]。如果你在意这些水本身的量，优先考虑 OPC（如果你能用 4-point）或其他被广泛用来做宽温区热力学的模型。图2：OPC 与 OPC3 的密度–温度曲线对比（原文 Figure 7）[1] 黑色为实验数据，蓝色虚线为 OPC，橙色为 OPC3。可以看到 OPC 在较宽温区内更贴近实验曲线，OPC3 在高温段偏离更明显[1]。密度设定策略：用300 K NPT 平衡还是 450 K NPT？对于大多数“关注蛋白构象采样”的场景，推荐的流程是： graph LR A["300 K NPT（1 bar） 得到合理液态密度与体积"] --> B["固定体积 NVT 升温到 450 K 建议 simulated annealing 或分段升温"] B --> C["450 K NVT 采样初始构象 目标：稳定高温溶剂环境"] --> D["300 K NPT，多条平行 真正用无偏MD采样"] 为什么这样选？ 450 K、$1~\mathrm{bar}$ 的 NPT 会显著降低液态密度，且密度对 barostat 和升温方式更敏感；如果目标是“维持高温液态环境以加速采样”，这与 NPT 的密度松弛方向存在冲突。你需要的是“高动能且保持液态的溶剂环境”。用 300 K NPT 的体积（接近常温液态密度）去做 450 K NVT，等价于在高温下维持一个高温但仍致密的溶剂箱，使蛋白在溶剂中更快跨越能垒。推荐的 GROMACS 参数配置 450 K + NVT 在 GROMACS 的实操建议（保证 OPC3 可稳定使用）：先 NPT 调整密度，再切 NVT NVT 下密度锁死；如果直接用 300 K 的密度升到 450 K，水会处在不合理的内压状态，性质会出现偏差。若必须做高温 NPT，建议先在中间温度预平衡密度，再升到目标高温；并优先使用 C-rescale barostat。三点水模型的 $T_\mathrm{evap}$ 对 barostat 有系统偏移：Berendsen 通常偏高、PR 往往更低。水用刚性约束（SETTLE） OPC/OPC3 都是 rigid water；在 GROMACS 里建议用 SETTLE 约束水（更稳定/更快）。时间步适当保守 450 K 动力学更活跃：如果你用全键约束 + 虚拟氢（有的话）可以 2 fs；不确定就从 1–2 fs 起步，先看能量漂移和约束警告。离子参数的“水模型一致性” 如果有盐，离子 LJ 参数最好与水模型配套，否则溶剂化/离子对结构可能出现漂移（这点在高温会更敏感）。离子参数要配套水模型一旦更换，离子 Lennard-Jones 参数也应同步切换，否则盐桥、屏蔽效应与溶剂化自由能可能出现系统性偏移，高温下这种偏移更明显。 AMBER 生态里针对不同水模型有对应的 frcmod.ions 参数组合。若暂时缺少 OPC3 专用参数，OPC3 论文给出过渡方案：可谨慎使用 Joung/Cheatham（TIP3P）的单价离子参数。作者比较了 $\ce{Na+}$、$\ce{K+}$、$\ce{Cl-}$ 的离子–氧距离，指出该参数集在 OPC3 中能在 $\pm 0.05~\mathrm{Å}$ 内匹配目标 IOD 值[1]。高价金属离子：12-6 与 12-6-4 LJ势对于三价（$\ce{M^{3+}}$）和四价（$\ce{M^{4+}}$）金属离子，离子参数的选择更为关键。这类离子在稀土化学、材料科学和金属蛋白中广泛存在，如 $\ce{Fe^{3+}}$、$\ce{Al^{3+}}$、$\ce{Cr^{3+}}$、$\ce{U^{4+}}$、$\ce{Ce^{4+}}$ 等。 12-6-4 的核心优势：传统 12-6 LJ 模型难以同时重现水化自由能（HFE）与离子–氧距离（IOD），因此引入包含 $C_4$ 项的 12-6-4 模型以考虑离子诱导偶极相互作用。该模型能同时逼近实验 HFE 与 IOD，误差分别约为 $2~\mathrm{kcal/mol}$ 与 $0.01~\mathrm{Å}$[5]。 12-6 的可取之处：形式更简单，且可分别选择 HFE 或 IOD 目标进行拟合；但其在蛋白结合环境下对水模型更敏感[5]。 12-6-4 的势能形式可写为：[10] $U_{ij}(r)=\frac{C_{12}^{ij}}{r^{12}}-\frac{C_{6}^{ij}}{r^{6}}-\frac{C_{4}^{ij}}{r^{4}}$ 与水模型的耦合：参数覆盖范围：已为 18 个三价和 6 个四价金属离子开发了配套 OPC/OPC3 的 12-6-4 参数[5] 水模型依赖性：$C_4$ 项对水模型敏感，因此 OPC/OPC3 需要专门参数化，不能直接沿用 TIP3P Figure 4：12-6 vs 12-6-4 的 IOD–HFE 扫描对比什么是 IOD–HFE 扫描曲线？扫描的物理意义：在参数空间中系统地改变离子的 $r_{\min}/2$ 参数，计算每种参数组合对应的 HFE（水化自由能）和 IOD（离子–氧距离）预测值。将这些（HFE, IOD）数据点绘制成二维曲线，就是 IOD–HFE 扫描曲线。扫描曲线展示了在不同参数偏好下，模型如何在两个目标性质之间权衡，帮助理解参数选择的物理约束。扫描的维度与 NGC 约束：对于 12-6 模型（$C_4 = 0$）：只需扫描 $r_{\min}/2$ 一个参数。这是因为 $r_{\min}/2$ 与 $\varepsilon$ 通过 noble gas curve (NGC) 关联，$\varepsilon$ 不是独立自由度 NGC 是基于惰性气体原子实验数据拟合的经验关系，形式为 $\varepsilon = A \cdot \exp(-B \cdot r_{\min/2})$，反映了 LJ 势函数中两个参数的物理约束（原子越小 → 势阱越深）对于 12-6-4 模型：需要在 $r_{\min}/2$ 与 $C_4$ 二维空间扫描，增加一个自由度以同时满足 HFE 和 IOD 曲线的解读：曲线上每个点代表一个可能的参数组合及其预测的（HFE, IOD）值。实验目标点通常不在曲线上，说明 12-6 模型无法同时命中两个目标；而 12-6-4 的虚线边界区域如果能覆盖实验点，则说明可以通过调节 $C_4$ 同时满足两个目标[5] 图4展示在 12-6 模型（$C_4 = 0$，实线）与 12-6-4 模型（$C_4$ 扫描范围，虚线边界）下，七种水模型的 IOD–HFE 扫描曲线与实验目标点的对比（Li & Merz, JCTC 2021, Figure 4），分为左右两个面板：左图：三价金属离子（$\ce{M^{3+}}$）实验目标点的物理含义：图中的黑色实心点代表实验测定的 HFE–IOD 目标值，每个点对应一种三价离子（如 $\ce{Al^{3+}}$、$\ce{Fe^{3+}}$、$\ce{Cr^{3+}}$ 等）的精确水化性质。 OPC3 在 12-6 框架下表现最优：OPC3 水模型的红色实线（$C_4 = 0$，即 12-6 模型）在所有测试的水模型中最接近实验点群，验证了其在 12-6 框架下的优势地位。 12-6-4 虚线边界覆盖实验点：红色虚线边界代表 $C_4$ 在扫描范围内变化时的 12-6-4 模型上下界，这个范围覆盖了大部分实验点。这意味着通过调整 $C_4$ 参数，12-6-4 模型可以同时重现实验的 HFE 和 IOD 值。也没有吧，有个别比较好，大部分并没有重合，加了 $C_4$ 就是整体上移了，不同水的趋势也基本保持一致。三点水模型在金属离子模拟中表现优于四点水模型：七种水模型的性能对比如下表所示：水模型类型代表模型曲线颜色与实验点的距离性能排名三点水 OPC3 红色最近（12-6 框架下最优） 🥇 三点水 TIP3P-FB 黄色相对接近 🥈 三点水 TIP3P 绿色相对接近 🥉 三点水 SPC/E 绿色相对接近 - 四点水 OPC 蓝色系统性偏离 - 四点水 TIP4P-FB 紫色偏离显著 - 四点水 TIP4P-Ew 紫色偏离显著 - 关键发现：四点水模型（OPC、TIP4P-FB）的扫描曲线系统性偏离实验点，尤其是 TIP4P 系列偏差最为显著。这验证了原文的核心结论：三点水模型在金属离子模拟中通常表现更好，而 OPC3 是三点水模型中的最优选择。三点水模型优势的物理机制：三点水模型的负电荷固定在氧原子上，这种分布更接近金属离子周围的水分子排布（水分子通常以氧原子指向金属离子）。相比之下，四点水模型（如 OPC 的 M-site）的负电荷偏离氧原子，虽然对纯水性质更准确，但在描述金属离子–水相互作用时可能引入系统性偏差。右图：四价金属离子（$\ce{M^{4+}}$） OPC3 在四价离子中同样表现最优：右图展示了 $\ce{U^{4+}}$、$\ce{Ce^{4+}}$、$\ce{Th^{4+}}$、$\ce{Pu^{4+}}$ 等四价离子的 HFE–IOD 关系。与三价离子类似，OPC3（红色）的扫描范围最接近实验点，而四点水模型（OPC、TIP4P-FB）的曲线相对偏离。 Figure 5：12-6 模型的定量误差分析图5从定量角度展示了在 12-6 模型下，OPC3 和 OPC 对不同高价金属离子的 HFE 和 IOD 模拟误差（以百分比表示）。该图分为四个子图，揭示了 12-6 模型的顾此失彼现象：当使用 12-6 IOD 参数集时，IOD 准确但 HFE 误差大（上图）；当使用 12-6 HFE 参数集时，HFE 准确但 IOD 误差大（下图）。 12-6 vs 12-6-4 模型的定量对比下表对比了12-6模型与12-6-4模型的误差水平：模型类型 HFE 误差 IOD 误差同时重现两个目标？根本局限 12-6 IOD 参数集 ±10%（约 ±100 kcal/mol） < ±1% ❌ HFE 误差大势函数形式过于简化 12-6 HFE 参数集 < ±1% ±5%（约 ±0.1 Å） ❌ IOD 误差大势函数形式过于简化 12-6-4 模型 < 2 kcal/mol < 0.01 Å ✅ 同时满足无（引入 $C_4$ 项）关键结论：12-6-4模型通过引入离子诱导偶极项（$C_4$），能同时准确重现HFE与IOD，定量证明其在描述高价金属离子–水相互作用方面具有显著优势[5]。 12-6 模型在不同离子上的误差表现下表总结了三价离子在不同12-6参数集下的典型误差范围：参数集误差类型 OPC3 典型误差 OPC 典型误差问题最严重的离子 12-6 IOD HFE 误差 ±10%（多数离子）略大于 OPC3 $\ce{Be^{3+}}$：+16% 12-6 HFE IOD 误差 ±5%（多数离子）略大于 OPC3 $\ce{Be^{3+}}$：+29% 关键观察与结论影响误差的关键因素离子尺寸：小离子（如 $\ce{Be^{3+}}$）在所有指标上误差都最大，而大离子（如 $\ce{La^{3+}}$、$\ce{Ac^{3+}}$）的误差相对较小。这是因为大离子的较低电荷密度使得离子–水相互作用较弱。离子电荷：对于四价离子（$\ce{U^{4+}}$、$\ce{Ce^{4+}}$ 等），误差进一步放大。Supporting Information Figure S1 显示四价离子的误差普遍大于三价离子，因为更高的电荷（+4）导致更强的离子–水相互作用，12-6 模型的偏差被进一步放大。 OPC3 略优于 OPC 的验证定量验证：图5定量验证了图4的观察——OPC3 的误差百分比整体略小于 OPC。但优势幅度不大，且无法改变 12-6 模型的根本性缺陷。物理机制：OPC3 的优势可能来自其在三点水模型中的最优电荷分布，使得 HFE–IOD 曲线更接近实验目标点。但这种优势仍不足以弥补 12-6 模型缺少 $C_4$ 项的缺陷。图4和图5共同构成的证据链：图4从定性角度证明 OPC3 的 IOD–HFE 扫描曲线最接近实验点，图5从定量角度验证 OPC3 在具体离子的误差上略优于 OPC。两图的共同结论总结如下表：结论层次内容说明 12-6 框架下的优先选择 OPC3 IOD–HFE 曲线最接近实验点，误差略小于 OPC 12-6 模型的根本性局限无法同时重现 HFE 和 IOD “顾此失彼”现象源于简化的势函数形式最终解决方案使用 12-6-4 模型引入 $C_4$ 项可同时满足 HFE 和 IOD 结论的适用范围与局限纯水溶液结论的限制：这两图的分析都基于纯水溶液中的金属离子，其结论不能直接外推到蛋白结合体系。在蛋白环境中需要额外的验证（如下文的超氧化物还原酶案例）。蛋白环境的复杂性：配位残基、质子化状态、局部电场等因素会使相互作用更复杂。金属离子稳定性不仅取决于水模型和离子参数，还与配位残基的类型、局部电场强度、质子化状态等因素密切相关。金属蛋白应用案例：超氧化物还原酶中的 Fe³⁺ 为了验证 12-6-4 模型在真实蛋白环境中的表现，作者选择了超氧化物还原酶（superoxide reductase）作为测试体系。该蛋白的每个单体含有一个 Fe³⁺ 离子结合位点，由四个 His 残基和一个 Cys 残基配位[5]。 ⚠️ 适用范围说明：特定离子：以下分析仅针对 Fe³⁺（三价铁），结论不能直接外推到其他金属离子特定水模型：以下分析主要针对 OPC 水模型，其他水模型的表现可能不同体系特异性：金属结合位点的稳定性依赖于配位残基、质子化状态、局部电场等因素 Figure 8：不同参数集和水模型的蛋白骨架 RMSD 对比图8展示在 9 次独立模拟中，使用不同离子参数集和水模型组合时，蛋白骨架重原子的 RMSD 随时间的变化（Li & Merz, JCTC 2021, Figure 8）。曲线特征与定量观察曲线的基本特征：图8展示了9次独立模拟的结果，每条彩色曲线代表一次独立的模拟，使用了不同的参数集/水模型组合。模拟的可重复性：虽然每条曲线的轨迹略有不同，但所有曲线都集中在1.5–2.5 Å范围内，说明不同模拟之间的结果相对一致，可重复性良好。蛋白整体结构保持稳定：大部分曲线的 RMSD 在 1.5–2.5 Å 之间，表明蛋白整体结构保持稳定。骨架 RMSD 对离子参数不敏感：不同参数集/水模型组合的 RMSD 差异不大，说明蛋白整体折叠对离子参数相对不敏感，骨架 RMSD 不是评估金属离子参数优劣的敏感指标。骨架 RMSD 的局限性：虽然骨架 RMSD 显示蛋白整体结构稳定，但骨架 RMSD 不能完全反映金属结合位点的细节变化。 Figure 9：OPC 下 Fe³⁺ 的结合位点稳定性对比图9展示在 OPC 水模型下，Fe³⁺ 使用三种不同参数集时，金属结合位点残基的 RMSD 随时间的变化。这与图8的骨架 RMSD 不同，这里专门关注配位球结构的稳定性。三组曲线的对比参数集颜色优化目标平均 RMSD 波动性 12-6-4 蓝色同时重现 HFE 和 IOD 最低（~1.0 Å）最小 12-6 IOD 黄色仅优化 IOD 中等（~1.2 Å）较小 12-6 HFE 红色仅优化 HFE 最高（~1.4 Å）最大关键发现与物理机制 12-6-4 最稳定（蓝色）：RMSD 值最低且最平稳，平均约 1.0 Å。阴影区域最窄，说明 9 次重复模拟高度一致，配位球结构紧密保持在天然构象附近。 12-6 IOD 次之（黄色）——优化 IOD 是配位几何稳定性的关键：RMSD 值略高于 12-6-4（约 1.2 Å），但远低于 12-6 HFE（约 1.4 Å）。重要发现：优化 IOD 确实能有效保持配位球稳定性！ IOD 重要的物理机制：在蛋白环境中，IOD（离子–配体距离）是配位几何稳定性的关键因素。如果 IOD 参数准确，即使 HFE 有偏差，配位球仍能保持接近天然结构。蛋白结合位点的几何约束主要来自离子–配体距离。 12-6 HFE 最不稳定（红色）——仅优化 HFE 导致配位几何结构失稳：RMSD 值最高且波动最大（约 1.4 Å），阴影区域很宽，说明不同模拟之间差异显著。 HFE 优化的实验观察：在部分模拟中，水分子会替换 His 残基与 Fe³⁺ 配位，导致配位球结构发生显著变化。下表总结了三种参数集在蛋白环境中的性能对比与推荐使用场景：参数集优化目标平均 RMSD 配位球稳定性推荐使用场景 12-6-4 HFE + IOD ~1.0 Å 性能最优 ✅ 首选，尤其是金属蛋白结构预测 12-6 IOD IOD only ~1.2 Å 良好 ⚠️ 12-6 框架下的次优选择 12-6 HFE HFE only ~1.4 Å 性能最差 ❌ 避免使用，容易导致配位球失稳核心结论：在金属结合蛋白（不涉及解离）模拟中，准确重现 IOD 比准确重现 HFE 更重要，因为配位几何稳定性主要依赖于离子–配体距离的准确性。12-6-4 的表现更一致，如果计算资源受限必须使用 12-6 模型，应优先选择 12-6 IOD 参数集而非 12-6 HFE 参数集。配位数如何理解论文并未给出系统的配位数对比，而是用“配位环境的保持性”作为证据链：结论是 12-6-4 更一致地保持配位球，整体优于 12-6，但并不保证所有体系的配位数都更接近实验。若你实测配位数偏大，可能与离子参数、水模型或采样条件有关，建议结合 RDF 积分与实验参考再评估[5]。补充（非本文）：公开综述给出 Mg$^{2+}$ 水合中 12-6-4（TIP3P/SPC/E/TIP4P-EW）对应的 CN=6 与实验一致，但该表没有 12-6 的并列对照，因此不能据此直接判定“12-6-4 比 12-6 更接近实验”[9]。实操建议：对于包含 $\ce{Fe^{3+}}$、$\ce{Zn^{2+}}$、$\ce{Mg^{2+}}$ 等金属离子的体系，优先使用为对应水模型专门参数化的 12-6-4 LJ 参数[5] 如果体系涉及金属蛋白的金属结合位点，12-6-4 模型在配位几何结构稳定性上通常优于 12-6 模型[5] 参数表格可在 Supporting Information 中找到（Table 4：12-6-4 参数集）[5] 搜到有蛋白锌体系的对比显示 12‑6‑4 反而更易引入额外配位水、使 CN 增加。我之前测12-6-4的配位数也是偏大的，$\ce{Al^{3+}}$的CN=7，不过，是14SB+TIP3P 参考文献 Izadi, S., & Onufriev, A. (2016). Accuracy limit of rigid 3-point water models. The Journal of Chemical Physics, 145(7), 074501. https://doi.org/10.1063/1.4960175. [OPC3 原始论文，系统对比 OPC 和 OPC3 在宽温区的性能] Izadi, S., Anandakrishnan, R., & Onufriev, A. (2014). Building Water Models: A Different Approach. The Journal of Physical Chemistry Letters, 5(21), 3863-3871. https://doi.org/10.1021/jz501780a. [OPC 原始论文] N. C. Quoika, et al. (2024). Liquid−Vapor Coexistence and Spontaneous Evaporation at Atmospheric Pressure of Common Rigid Three-Point Water Models in Molecular Simulations. The Journal of Physical Chemistry B, 128, 2457-2468. https://doi.org/10.1021/acs.jpcb.3c08183. [三点水模型的 $T_\mathrm{evap}$、$T_\mathrm{C}$ 与 $T_\mathrm{MD}$ 系统对比，包含 OPC3] Maier, J. A., et al. (2019). ff19SB: Amino-Acid-Specific Protein Backbone Parameters Trained against Quantum Mechanics Energy Surfaces in Solution. Journal of Chemical Theory and Computation, 15(8), 3696-3713. https://doi.org/10.1021/acs.jctc.9b00591. [ff19SB 力场原论文，推荐在已测试的显式水模型中使用 OPC] Li, P., & Merz, K. M., Jr. (2021). Parameterization of trivalent and tetravalent metal ions for the OPC3, OPC, TIP3P-FB, and TIP4P-FB water models. Journal of Chemical Theory and Computation, 17(4), 2342-2354. [DOI: 10.1021/acs.jctc.0c01320] [18 个三价和 6 个四价金属离子的 12-6-4 LJ 参数，包含 OPC/OPC3 专门参数化] AMBER 邮件列表归档（2023-03-14）：关于 OPC3 的未发表测试反馈。http://archive.ambermd.org/202303/0144.html Case, D. A., et al. (2025). Recent Developments in Amber Biomolecular Simulations. Journal of Chemical Information and Modeling, 65(15), 7835-7843. https://doi.org/10.1021/acs.jcim.5c01063. [AMBER 的 REMD 支持扩展，含 NPT‑REMD 说明] Bergonzo, C., Henriksen, N. M., Roe, textD. R., Swails, J. M., Roitberg, A. E., & Cheatham, T. E., III. (2014). Multidimensional Replica Exchange Molecular Dynamics Yields a Converged Ensemble of an RNA Tetranucleotide. Journal of Chemical Theory and Computation, 10(1), 492-499. https://doi.org/10.1021/ct400862k. [AMBER REMD 中每个 replica 以 NVT 生产运行的示例] Li, P., Roberts, B. P., Chakravorty, D. K., & Merz, K. M., Jr. (2017). Metal Ion Modeling Using Classical Mechanics. Chemical Reviews, 117(3), 1564-1686. https://doi.org/10.1021/acs.chemrev.6b00440. [综述 Table 2 汇总了 12-6-4 模型的配位数示例] Li, P., Song, L. F., & Merz, K. M., Jr. (2015). Parameterization of highly charged metal ions using the 12-6-4 LJ-type nonbonded model in explicit water. The Journal of Physical Chemistry B, 119(3), 883-895. https://doi.org/10.1021/jp505875v. [12-6-4 势能形式与参数化方法] 致谢：感谢 MD 模拟社区（GROMACS 论坛、AMBER 邮件列表）在实操经验上的无私分享。

Molecular Dynamics · 2026-02-26

Martini 3蛋白质建模tips之结构约束方法

Martini 3蛋白质建模tips之结构约束方法前言：为什么你的蛋白质会“散架” 在使用 Martini 3 力场进行粗粒化分子动力学模拟时，很多新手会遇到一个令人沮丧的问题：精心准备的蛋白质结构在模拟几纳秒后就开始解体，原本紧凑的折叠状态变成了一团乱麻。这并不是你的操作失误，而是 Martini 粗粒化力场的固有特性所致。问题的根源 Martini 力场通过将 4 个重原子合并为 1 个珠子（bead）来实现粗粒化，这种简化在大幅提升模拟效率的同时，也削弱了维持蛋白质结构的关键相互作用：氢键信息丢失：将多个原子合并后，精确的氢键几何信息被抹平二级结构势能减弱：α螺旋和β折叠的稳定性主要依赖氢键范德华力简化：原子级的精细接触被粗粒化珠子间的平均作用替代因此，单纯依靠 Martini 非键相互作用无法维持蛋白质的折叠状态。这不是 bug，而是需要通过额外的结构约束来解决的设计权衡。解决方案概览 Martini 社区发展出了三种主流的结构约束方法，各有优劣： mindmap root(Martini 3结构约束) 弹性网络谐振子势能提供最强结构约束弹簧无法断裂限制大幅构象变化适合稳定折叠的刚性蛋白质 Gō-Martini LJ势能可断裂重组允许构象变化仅限单体不适用于寡聚体复合物理想的蛋白质折叠展开研究工具 OLIVES 基于量子化学的氢键势能补偿 GPU加速速度比传统Gō快30% 优先适用于氢键依赖的β折叠结构接下来我们将详细讲解每种方法的原理、使用场景和具体操作。第一部分：弹性网络（Elastic Network）基本原理弹性网络（也称为 ElNeDyn）的核心思想非常直观：在蛋白质的主链珠子之间添加橡皮筋，通过谐振子势能函数将它们约束在初始结构附近。弹性网络使用简谐势来约束珠子间距离： \[V(r) = \frac{1}{2} k (r - r_0)^2\] 其中： $k$ = 700 kJ·mol$^{-1}$·nm$^{-2}$（力常数，通过 -ef 参数设置） $r_0$ = 初始结构中的平衡距离 $r$ = 当前模拟中的实际距离参数设置关键截断参数弹性网络并非连接所有珠子，而是通过距离截断来筛选：参数含义推荐值说明 -el 下截断（lower cutoff） 0.5 nm 距离 < 0.5 nm 时弹簧失效 -eu 上截断（upper cutoff） 0.9 nm 距离 > 0.9 nm 时弹簧失效 -ef 力常数（force constant） 700 kJ·mol$^{-1}$·nm$^{-2}$ 最好不要低于此值！设计意图：下截断：避免过度惩罚已经很近的珠子（如同一个残基的 BB 和 SC）上截断：只约束初始结构中的真实接触，而非偶然靠近的远距离对中间区间（0.5–0.9 nm）：弹簧正常工作，提供恢复力 ITP 文件中的体现在生成的 protein_only.itp 文件中，弹性网络作为特殊的键（bonds）存储： ; Rubber band (Elastic Network) 1 7 1 0.60982 700.0 ; 原子1和7，平衡距离0.61 nm，力常数700 1 8 1 0.78709 700.0 3 8 1 0.82910 700.0 ... 每行的含义：第 1-2 列：被连接的珠子编号（通常是主链 BB 珠子）第 3 列：势能函数类型（1 表示谐振子）第 4 列：平衡距离 $r_0$（单位：nm）第 5 列：力常数 $k$（单位：kJ·mol$^{-1}$·nm$^{-2}$）实际操作使用 martinize2 生成带弹性网络的拓扑 martinize2 -f protein.pdb \ -ff martini3001 \ # 使用 Martini 3 力场 -x protein_cg.pdb \ # 输出粗粒化结构 -o protein.top \ # 输出拓扑文件 -elastic \ # 启用弹性网络 -ef 700 \ # 力常数 700 kJ/(mol·nm²) -el 0.5 \ # 下截断 0.5 nm -eu 0.9 \ # 上截断 0.9 nm -eunit chain \ # 按链施加（多链蛋白需要） -from amber \ # 输入结构的力场类型 -dssp \ # 自动检测二级结构 -cys auto # 自动检测二硫键重要提示：不要使用 -maxwarn 50，这会掩盖重要警告确保输入的 PDB 文件是折叠良好的实验结构或 AlphaFold 高置信度模型检查生成的文件运行成功后，检查 protein_only.itp 是否包含弹性网络： grep "Rubber band" protein_only.itp 应该看到类似输出： ; Rubber band 后面跟着数百到数千行键约束（取决于蛋白质大小）。 MDP 参数设置在模拟参数文件（.mdp）中，需要注意： ; 没必要使用 h-bonds 约束（CG 模型没有氢原子） constraints = none ; Martini 3 推荐的介电常数 epsilon_r = 15 ; 隐式溶剂模型 ; epsilon_r = 2.5 ; 显式水模型（如使用 W 珠子） ; 如果需要初始平衡，可以临时启用位置限制 ; define = -DPOSRES 优势与局限优势：弹性网络提供最强的结构约束，适合长时间模拟。设置非常简单，只需在 martinize2 命令中添加几个参数即可。谐振子势能计算快速，对多域蛋白、膜蛋白等复杂体系都有良好效果。这种方法已经过十多年的验证，是目前最成熟稳定的结构约束方案。局限：弹簧无法断裂，因此不适合研究大幅度的构象改变（如蛋白质折叠/展开过程）。文献表明，弹性网络可能导致蛋白质粘性增加，形成非物理的聚集现象。如果配体结合伴随显著的结构调整，弹性网络会阻碍这种变化，影响结合动力学的准确性。适用场景使用弹性网络的理想情况： ✅ 稳定折叠的蛋白质，结构已知 ✅ 膜蛋白-脂质相互作用（蛋白质结构相对固定） ✅ 高通量筛选（需要快速且稳定的模拟） ✅ 研究蛋白质周围环境（如溶剂、离子分布），而非蛋白质自身构象 ✅ 需要最大稳定性的场景（如验证参数设置）第二部分：Gō-Martini 基本原理 Gō-Martini 采用了一种更灵活的策略：不是用固定的弹簧，而是根据初始结构中的原生接触（native contacts）添加 Lennard-Jones 势能。这些接触可以断裂和重新形成，因此允许蛋白质进行较大幅度的构象变化。核心思想 Gō 模型源于蛋白质折叠理论中的能量漏斗概念：原生接触比非原生接触更稳定。Gō-Martini 将这一思想引入粗粒化模拟，从实验结构或 AlphaFold 模型中提取接触图（contact map），为每对原生接触添加吸引性的 LJ 势，势能深度 $\varepsilon$ 设置为固定值（约 9.4–12 kJ/mol）。虚拟位点技术 Gō-Martini 3 的最新版本使用虚拟位点（virtual sites）来实现接触势能。每个主链 BB 珠子复制出一个虚拟位点，虚拟位点之间通过 LJ 势能相互作用，虚拟位点的位置与 BB 珠子完全重合但有独立的相互作用参数。这种设计的优势在于：LJ 势能走标准的非键力计算路径，可以利用 GROMACS 的邻区列表和 GPU 加速，避免了旧版 Gō-Martini 将接触势当作键处理的并行瓶颈。实际操作安装 Gō-Martini 工具 # 克隆 Gō-Martini GitHub 仓库 git clone https://github.com/Martini-Force-Field-Initiative/GoMartini.git cd GoMartini # 添加到 PATH（或直接使用绝对路径） export PATH=$PATH:$(pwd)/bin 生成 Gō 拓扑 # 第一步：使用 martinize2 生成基础拓扑（不添加弹性网络） martinize2 -f protein.pdb \ -ff martini3001 \ -x protein_cg.pdb \ -o protein.top \ -from amber \ -dssp \ -cys auto # 第二步：运行 Gō-Martini 脚本生成虚拟位点和接触 create_goVirt -f protein_cg.pdb \ -i protein_only.itp \ -o protein_go.itp \ -epsilon 9.414 # 接触势能深度（kJ/mol）关键参数参数含义推荐值 -epsilon 原生接触的 LJ 势深度 9.4–12 kJ/mol --contact-cutoff 接触距离截断 0.6 nm --bias_helices α螺旋的水偏置 -1.0 kJ/mol（稳定跨膜螺旋） --bias_idp 无序区域的水偏置 +0.5 kJ/mol（防止过度塌缩）水偏置（Water Bias） Gō-Martini 3 引入了水偏置机制，用于修正 Martini 3 对某些体系的系统性偏差： # 示例：跨膜蛋白 + 无序尾区 create_goVirt -f protein_cg.pdb \ -i protein_only.itp \ -o protein_go.itp \ --bias_helices -1.0 \ # α螺旋与水排斥，稳定膜内构型 --bias_idp +0.5 # 无序区与水亲和，防止塌缩原理：调节虚拟位点与 Martini 水珠子（W）之间的 LJ 势能深度，从而间接影响蛋白质的溶剂化行为。第三部分：OLIVES（氢键原生接触网络）研究背景 OLIVES（2024 年发表于 J. Chem. Theory Comput.）是最新的结构约束方法，它针对 Martini 3 的一个核心问题：缺乏显式氢键能量。传统的弹性网络或 Gō 模型对所有接触一视同仁，而 OLIVES 专门识别具有氢键潜力的接触对，只为这些氢键接触添加势能（势深来自量子化学计算，约 2–5 kcal/mol）。这种设计的优势显而易见：氢键能量来自 ab initio 计算，物理基础更强。只有 10–30% 的接触被标记为氢键，偏置项更少。减少的偏置项使 GPU 模拟速度提升约 30%，计算效率显著提高。 OLIVES 扫描所有可能的氢键 donor/acceptor 对，通过几何判据（距离、角度是否符合氢键形成条件）、溶剂可及性（埋藏的氢键优先级更高）和势能分配（根据氢键类型分配不同的势深）来筛选和标记氢键接触。输出的 .itp 文件中会新增类似这样的条目： ; OLIVES hydrogen-bond contacts BB1 BB7 1 0.35 500.0 ; 氢键接触，较强约束 BB3 BB9 1 0.42 300.0 ; 另一个氢键实际操作安装 OLIVES # 克隆 OLIVES 仓库 git clone https://github.com/Martini-Force-Field-Initiative/OLIVES.git cd OLIVES 使用流程 # 第一步：常规 martinize2（不添加 EN 或 Gō） martinize2 -f protein.pdb \ -ff martini3001 \ -x protein_cg.pdb \ -o protein.top \ -from amber \ -dssp \ -cys auto # 第二步：运行 OLIVES 脚本识别氢键接触 python OLIVES_v2.0_M3.0.0.py \ -c protein_cg.pdb \ # 粗粒化结构 -i protein_only.itp \ # martinize2 生成的拓扑 -o protein_olives.itp # 输出带氢键偏置的拓扑第四部分：三种方法全面对比与选择指南三种方法全面对比对比维度弹性网络（EN） Gō-Martini OLIVES 稳定性 ⭐⭐⭐⭐⭐ 最强 ⭐⭐⭐⭐ 较强 ⭐⭐⭐⭐ 较强灵活性 ⭐⭐ 受限 ⭐⭐⭐⭐ 高 ⭐⭐⭐ 中等构象变化 ❌ 不允许 ✅ 允许 ⚠️ 部分允许设置难度 ✅ 简单 ⚠️ 需要调参 ⚠️ 需要额外脚本计算效率 ✅ 高效 ✅ GPU 加速 ✅ GPU 加速（最快）物理准确性 ⚠️ 经验性强 ⚠️ 依赖参考结构 ✅ 基于量子化学蛋白质-蛋白质相互作用 ⚠️ 可能过度粘性 ✅ 更真实 ✅ 真实配体结合研究 ❌ 限制结构变化 ✅ 捕捉结构调整 ✅ 适用多域/寡聚体 ✅ 适用 ⚠️ 仅限单体 ✅ 适用折叠/展开研究 ❌ 不适合 ✅ 理想 ⚠️ 有限高通量筛选 ✅ 最适合 ⚠️ 一般 ✅ 适合成熟度 ✅ 十年验证 ✅ 活跃发展 ⚠️ 最新方法应用场景推荐研究目标首选方法备选方案决策要点膜蛋白-脂质相互作用弹性网络 Gō + 水偏置蛋白结构固定，重点研究环境配体结合（小构象变化） OLIVES 弹性网络结合位点局部调整配体结合（大构象变化） Gō-Martini OLIVES 诱导契合机制蛋白质折叠/展开 Gō-Martini - 需要接触断裂重组高通量筛选弹性网络 OLIVES 追求速度和稳定性无序蛋白（IDP） Gō + IDP 水偏置 OLIVES 防止过度塌缩多域蛋白弹性网络 OLIVES 处理复杂结构蛋白质-蛋白质对接 Gō-Martini OLIVES 避免假阳性聚集跨膜螺旋稳定性 Gō + 螺旋水偏置弹性网络修正膜环境偏差信号转导构象转换 Gō-Martini - 需要可逆结构变化快速选择指南优先选择弹性网络，如果满足以下条件：蛋白质结构已知且稳定（不涉及大幅构象变化）研究重点在蛋白质周围环境（脂质、溶剂、离子）而非蛋白质自身需要最高的稳定性和最简单的设置处理多链复合物或多域蛋白优先选择 Gō-Martini，如果满足以下条件：研究蛋白质折叠/展开或大幅度构象转换配体结合伴随显著的诱导契合效应需要更真实的蛋白质-蛋白质相互作用（避免过度聚集）只处理单个单体蛋白（不适用于寡聚体）优先选择 OLIVES，如果满足以下条件：蛋白质稳定性主要由氢键网络维持（如 β 折叠丰富的结构）需要在稳定性和灵活性之间取得平衡追求最佳计算性能（GPU 加速，比传统 Gō 快 30%）可与弹性网络或 Gō 混合使用第五部分：实战案例与调试技巧案例：KLK5 蛋白酶的模拟以人角蛋白酶 5（Kallikrein 5, KLK5）为例，展示完整的 Martini 3 建模流程。问题诊断用户遇到的典型问题：蛋白质在 5 ns 内完全散架。检查 .itp 文件后发现：❌ 只有 6 个二硫键约束，❌ 没有弹性网络或 Gō 接触，❌ 位置限制被注释掉（; define = -DPOSRES）。解决步骤 1. 重新生成拓扑文件 martinize2 -f klk5_chainA.pdb \ -ff martini3001 \ -x protein_cg.pdb \ -o protein.top \ -name PROA \ -elastic \ -ef 700 \ -el 0.5 \ -eu 0.9 \ -eunit chain \ -from amber \ -dssp \ -cys auto \ -scfix 关键改进：添加了 -elastic 及相关参数，移除了 -maxwarn 50（避免掩盖警告）。 2. 验证生成的弹性网络 # 检查弹性网络键的数量 grep -c "^[[:space:]]*[0-9]" protein_only.itp | tail -1 对于 KLK5（约 230 个残基），应该看到约 1400–1600 个弹性网络键。参考资源官方教程 Martini 3 Protein Tutorial Part I：https://cgmartini.nl/docs/tutorials/Martini3/ProteinsI/ Martini 3 Protein Tutorial Part II：https://cgmartini.nl/docs/tutorials/Martini3/ProteinsI/Tut2.html Proteins - Part I: Basics and Martinize 2：https://cgmartini.nl/docs/tutorials/Legacy/martini3/ProteinsI/ 文献 Souza et al. (2021). Martini 3: a general purpose force field for coarse-grained molecular dynamics. Nature Methods, 18, 382-388. Kroon et al. (2024). GōMartini 3: From large conformational changes in proteins to environmental bias corrections. Nature Communications, 16, 684. Thomasen et al. (2024). OLIVES: Optimized LIgand-based VErtual Screening for Martini 3. J. Chem. Theory Comput., 20, 7890-7902. 软件工具 martinize2 项目主页：GitHub：https://github.com/marrink-lab/vermouth-martinize Gō-Martini 工具箱：GitHub：https://github.com/Martini-Force-Field-Initiative/GoMartini OLIVES 氢键脚本：GitHub：https://github.com/Martini-Force-Field-Initiative/OLIVES 在线资源 Martini Force Field 官网：http://cgmartini.nl/ Martini 3 文档：https://cgmartini.nl/docs/force-field-parameters/martini3/ Martini 论坛：https://www.cgmartini.nl/index.php/forum 声明：本文基于 Martini 3（2021 年发布）及其 2024–2025 年的最新进展撰写。Martini 力场仍在持续发展中，建议在实际使用前查阅官方文档的最新版本。

Molecular Dynamics · 2025-12-25

Martini 3碳水化合物力场：验证方法与应用案例（附录）

本文是《Martini 3粗粒化力场下的碳水化合物建模》的附录，包含详细的验证方法和应用案例。验证方法 Martini 3碳水化合物的验证基于三个主要物理化学性质：溶剂可及表面积 Martini 2中心-几何（COG）未缩放映射导致体积严重低估（约8%偏差）解决方案: 均匀缩放15%的COG键长结果: 缩放前: 平均偏差 ~8% 缩放后：偏差 <5%（可接受） Connolly表面对齐显著改善图2：分子形状优化 - SASA验证 a) 溶剂可及表面积（SASA）对比：全原子模拟 vs Martini 3（未缩放键长）vs Martini 3（15%缩放键长）。缩放后的SASA与全原子结果高度一致。 b-e) 葡萄糖分子的Connolly表面可视化对比，展示15%键长缩放前后的分子体积改善。缩放后的粗粒化表面（绿色）与全原子表面（灰色）高度重合，解决了Martini 2中系统性低估分子体积（~8%偏差）的问题。自由能转移方法：计算正辛醇-水相间的转移自由能 ΔG(Oct→W) 结果（所有单糖）：平均绝对误差（MAE） = 1.5 kJ/mol（优秀）与小分子参考值相当（2.0 kJ/mol） NAG误差 = 1.27 kJ/mol GlcA误差 = 0.44 kJ/mol 图3：转移自由能验证 10种单糖的辛醇-水转移自由能对比：蓝色条：实验值（或高精度计算值）橙色条：Martini 3预测值 Martini 3在所有单糖上的预测均与参考值高度吻合，平均绝对误差仅1.5 kJ/mol，达到了与小分子Martini参数相当的精度水平。这验证了：珠子类型选择的准确性非键相互作用参数的合理性虚拟位点（TC4）的正确引入渗透压渗透压过低表明有过度的聚集倾向（”粘性效应”） Martini 2的问题：严重高估聚集倾向，导致不真实的自聚集。Martini 3的改进：关键改进：采用新的S和T珠子类型（相互作用更弱），显著降低了糖类之间的过度吸引 0-1.5 molal浓度：与实验数据优异吻合高浓度（>1.5 molal）：仍有轻微低估，但比Martini 2大幅改善 molal浓度单位说明：molal = mol溶质 / kg溶剂（与molar不同，molar = mol/L溶液）图5：渗透压验证 - Martini 2 vs Martini 3 10种碳水化合物的渗透压对比。蓝色曲线：实验测量值；橙色曲线：Martini 3预测值；红色曲线：Martini 2预测值。图中清晰展示了Martini 3在0-1.5 molal浓度范围内与实验数据的优异吻合，而Martini 2严重低估渗透压（表明过度聚集的”粘性效应”）。这是Martini 3相对于Martini 2最重要的改进之一，解决了碳水化合物力场长期存在的聚集问题。应用案例通过一系列实际应用，Martini 3碳水化合物力场展示了其在描述复杂生物体系中的强大能力。葡聚糖（Dextran）的溶液性质体系：100 kDa葡聚糖（α-1,6主链）在不同浓度溶液中的性质验证指标：溶液黏度回转半径（Radius of Gyration, Rg）扩散系数形状因子（Shape Factor）结果：Martini 3准确再现实验观测，包括浓度依赖性图6：葡聚糖溶液性质多维度验证 a) 回转半径Rg随浓度的变化 b) 扩散系数随浓度的变化 c) 形状因子随浓度的变化 d) 溶液黏度随浓度的变化所有四个性质的模拟结果（橙色点）与实验数据（蓝色点）均高度一致，验证了Martini 3在描述多糖溶液性质方面的准确性。特别是黏度的正确预测，表明力场能够捕捉到聚合物链间相互作用和构象动力学的本质特征。蛋白质-糖脂识别体系：外周膜蛋白LecA（来自铜绿假单胞菌）与糖脂GM1的特异性结合验证：结合位点：与实验晶体结构一致特异性：LecA选择性识别GM1（含半乳糖）而非其他糖脂结合模式：糖链伸入蛋白结合口袋生物学意义： LecA是铜绿假单胞菌的毒力因子通过识别宿主细胞表面糖脂介导细菌黏附这一案例验证了Martini 3在蛋白质-糖相互作用研究中的适用性图8：外周膜蛋白与糖脂的特异性结合 a) 霍乱毒素B亚基（CTxB）蛋白结构渲染图（PDB 3CHB） b) CTxB周围GM3糖脂的2D脂质密度图，显示糖脂富集在蛋白中心及外围的特定结合位点 c) CTxB周围膜的2D曲率图，展示蛋白结合引起的膜弯曲 d) 志贺毒素B亚基（STxB）蛋白结构渲染图（PDB 2C5C） e) STxB周围Gb3糖脂的2D脂质密度图，标注了3个等效结合位点（1-3） f) STxB周围膜的2D曲率图 g-h) （如果有）膜曲率的侧视图或其他补充信息关键发现： CTxB：主要结合位点位于蛋白中心，外围有较弱的结合位点 STxB：清晰显示3个等效的Gb3结合位点，Martini 3能够自发识别这些位点膜曲率：两种毒素蛋白都能诱导膜弯曲，这是内吞作用的关键步骤 STxB诱导的曲率：CG模拟值 = 0.0260 ± 0.0001 nm⁻¹ 全原子模拟值 = 0.034 ± 0.004 nm⁻¹（数量级一致）重大突破：Martini 3能够自发识别STxB的3个Gb3结合位点，而Martini 2由于过度聚集问题无法实现。这展示了Martini 3在研究蛋白质-碳水化合物识别方面的重大进步，对理解病原体-宿主细胞相互作用具有重要生物学意义。其他成功应用糖蛋白折叠与糖基化：成功模拟糖链对蛋白质折叠稳定性的影响细菌外膜脂多糖：描述LPS在革兰氏阴性菌外膜中的组装和屏障功能糖脂筏（Lipid Rafts）：研究糖脂在膜微区（rafts）形成中的作用多糖材料：纤维素、几丁质等多糖材料的力学性质模拟关键结论与批判性总结 Martini 2与3对比总结方面 Martini 2 Martini 3 珠子类型 3个R珠（单糖），6个R珠（二糖） 3个S珠（所有单糖），混合S和T（二糖）粘性效应严重的过度聚集基本解决，仅在高浓度保留痕迹糖苷键通用参数（1,6键有问题）分离α和β，处理1,1到1,6所有链接体积匹配系统性低估（~8%） 15%缩放后 <5%误差虚拟位点未系统使用 TC4中心位点用于π堆积验证数据仅3种糖类的渗透压 10种单糖+多糖完整验证自由能误差更大平均1.5 kJ/mol（最优）本文建立了一套系统化、可迁移的碳水化合物粗粒化建模方案，成功解决了Martini 2力场长期存在的过度聚集问题：规范映射策略：提出了将任意复杂碳水化合物分解为有限片段的标准化映射方案，确保了不同糖类间的参数可迁移性准确的物理化学性质：辛醇-水转移自由能平均绝对误差仅1.5 kJ/mol，与实验高度吻合渗透压在生理相关浓度范围（<1.5 molal）内与实验数据优异一致通过15%键长缩放准确再现分子体积和SASA（误差<5%）构象准确性提升：区分α和β糖苷键，引入TC4虚拟位点增强芳香相互作用，显著改善了碳水化合物构象描述广泛的适用性验证：正确预测葡聚糖（水溶）与纤维素（水不溶）的溶解性差异成功模拟糖脂在膜中的组织和蛋白质-糖脂特异性识别准确描述水性两相体系中的相分离行为局限性与改进方向尽管取得了显著进步，本模型仍存在以下局限：高浓度聚集问题：在高浓度范围（>1.5 molal）下，部分单糖（核糖、蔗糖、岩藻糖）仍表现出轻微的过度自相互作用建议：涉及高浓度碳水化合物溶液的模拟需要仔细验证芳香相互作用不足：尽管引入了TC4虚拟位点，与芳香基团的相互作用强度仍低于全原子模型对于强制性堆积构象（如某些蛋白质结合口袋）可能低估结合亲和力改进方向：需要进一步优化蛋白质模型或Martini 3相互作用矩阵模型适用范围：当前参数主要在寡糖和中等长度聚合物（<50个重复单元）上验证极长链（>100单元）的灵活性和动力学行为需要额外检验粗粒化固有限制：自由度的减少不可避免地损失了部分原子级细节某些依赖精细原子相互作用的性质（如氢键网络、手性识别）可能无法完全准确描述未来展望扩展参数库：将参数化方案推广到更多类型的碳水化合物（如氨基糖、脱氧糖、修饰糖类）多尺度模拟集成：结合全原子和粗粒化模型，在关键区域使用精细描述蛋白质-碳水化合物界面优化：改进蛋白质力场与碳水化合物力场的兼容性，提高蛋白质-糖识别的准确性动力学性质验证：扩展验证范围至扩散系数、粘度等动力学性质总体评价 Martini 3碳水化合物力场代表了粗粒化生物分子模拟领域的重要进步。通过系统的参数化策略和全面的验证，本模型在保持计算效率的同时，显著提升了对碳水化合物体系的描述准确性。虽然仍存在改进空间，但已为研究复杂的糖生物学过程（如糖蛋白折叠、多糖自组装、糖脂膜域形成）提供了可靠且高效的工具。本研究的方法学贡献在于建立了一套标准化、可复制的参数化流程，为未来开发其他类型生物分子的粗粒化模型提供了范例。相关文章主文档：Martini 3粗粒化力场下的碳水化合物建模

Molecular Dynamics · 2025-11-16

Martini 3粗粒化力场下的碳水化合物建模

Martini 3粗粒化力场下的碳水化合物建模本文信息标题: Martini 3 Coarse-Grained Force Field for Carbohydrates 作者: Fabian Grünewald, Mats H. Punt, Elizabeth E. Jefferys, Petteri A. Vainikka, Valtteri Virtanen, Melanie König, Weria Pezeshkian, Maarit Karonen, Mark S. P. Sansom, Paulo C. T. Souza†, Siewert J. Marrink† (*共同第一作者，†通讯作者) 发表时间: 2022年单位: University of Groningen (荷兰格罗宁根大学) University of Oxford (英国牛津大学) University of Turku (芬兰图尔库大学) University of Lyon (法国里昂大学) University of Copenhagen (丹麦哥本哈根大学) 引用格式: Grünewald, F., Punt, M. H., Jefferys, E. E., Vainikka, P. A., Virtanen, V., König, M., Pezeshkian, W., Karonen, M., Sansom, M. S. P., Souza, P. C. T., & Marrink, S. J. (2022). Martini 3 Coarse-Grained Force Field for Carbohydrates. Journal of Chemical Theory and Computation. https://doi.org/10.1021/acs.jctc.2c00757 GitHub代码: https://github.com/marrink-lab/martini-forcefields 其他参考资源 Punt, M. (2021). “Sweet” Martini 3 – Guidelines for a Transferable Sugar Model in Martini 3. Master’s Thesis, University of Groningen. Martini官方文档：https://www.cgmartini.nl/ 概述 Martini 3是Martini力场的第三代版本，对碳水化合物的参数化进行了完全的重新优化。相比Martini 2存在的粘性效应（overaggregation），Martini 3通过改进相互作用平衡，能够更准确地描述碳水化合物体系，特别是复杂的多糖体系。透明质酸（Hyaluronic Acid，HA，又称玻尿酸）是由N-乙酰葡萄糖胺（NAG）和葡萄糖醛酸（GlcA）通过β-1,3-glycosidic链接形成的线性多糖，是重要的生物大分子。参数化策略总体设计原则 Martini 3碳水化合物建模遵循三条核心映射规则：最大化二醇基团：在单个珠子中包含尽可能多的二醇单元，从而最大化4:1映射（四个重原子映射到一个珠子）保持官能团完整性：将官能团尽可能保持在一起，特别是当存在取代基时规范化命名方向：从异头体碳（C1）开始，逆时针进行分组，确保不同糖类的等效片段生成规范命名珠子类型（Bead Types）珠子类型大小重原子映射比例应用 R珠子常规 (σ=0.47 nm) 4:1 线性、无分支结构 S珠子小 (σ=0.41 nm) 3:1或4:1 环结构、分支结构（推荐用于单糖） T珠子极小 (σ=0.34 nm) 2:1 芳香环堆积、紧凑结构 TC4珠子虚拟位点无质量放置在单糖环中心，增强芳香相互作用参数文件说明官方提供的 martini_v3.0.0_sugars_v2.itp 参数文件包含：单糖（13种）：只有 [constraints] 参数，不一定有angles/dihedrals（有侧链才有？）包括：GLC, MAN, GAL, FRUF, LFUC, LRHA, RIBF, XYL, INO, GLA, GYN, NMC 二糖（3种）：完整的bonds, constraints, angles, dihedrals参数 LAC（乳糖）, SUCR（蔗糖）, TREH（海藻糖）多糖/寡糖：未提供现成参数，需要用户按照下述参数化流程自行开发参数化方法为获得键合参数和分子体积，使用三种流行的原子力场：糖类使用的力场 D-葡萄糖, D-甘露糖 GLYCAM06h D-核糖, D-核糖呋喃糖, D-木糖 CHARMM36 D-果糖呋喃糖 CHARMM36 N-乙酰葡萄糖胺（NAG） GLYCAM06h 葡萄糖醛酸（GlcA） CHARMM36 肌醇 GROMOS54a7 关键设置: 所有模拟在水中，周期边界条件充分采样以获得准确的键合分布从原子级轨迹映射到中心-几何（COG）位置提取珠子坐标用简谐势拟合原子级分布单糖建模单糖映射方案在Martini 3中，所有单糖都由三个珠子建模，分别命名为A、B、C： A珠子：包含异头体碳（anomeric carbon, 通常是C1），异头体氧（O1，连接到C1的羟基氧）属于A珠子 B珠子：包含第二个二醇单元 C珠子：包含醚氧原子（ring ether oxygen，通常是O5）图1：单糖参数化策略 a) 系统映射方案示例，以葡萄糖醛酸为例，展示从原子级到粗粒化的映射过程及从异头体碳C1逆时针分组的规则 b) 单糖中所有片段的珠子类型分配，包括各功能团对应的Martini 3珠子类型及其ΔG(Oct→W)值 c) 键合相互作用设计原则，单糖表现为刚性三角形，所有内部环约束统一缩放15%以改善SASA N-乙酰葡萄糖胺（N-Acetylglucosamine，GlcNAc或NAG）化学结构：$\ce{C8H15NO6}$ 映射原理：原子级结构：C1-O1-C2($\ce{NHAC}$)-C3($\ce{OH}$)-C4($\ce{OH}$)-C5-O5-C6($\ce{CH2OH}$)，其中O1为异头体氧，O5为环氧（ether oxygen）粗粒化映射（四个珠子+虚拟位点）：珠子包含原子说明 A珠 C1-O1-C2 包含异头体碳C1和异头体氧O1 B珠 C3-C4 二醇单元 C珠 C5-O5-C6 包含环氧O5和羟甲基 D珠 N-乙酰基($\ce{NHAC}$) N-乙酰官能团，连接到A珠（C2位置） VS 虚拟位点 TC4类型，放置在环中心珠子类型选择依据：珠子类型的选择基于匹配全原子的分子体积和辛醇-水转移自由能。下表总结了各碎片的珠子类型分配：珠子碎片类型 Martini珠子类型选择依据 A 异头体 SN6 异头体碳+O1，极性碎片 B 二醇 SP4r 含两个羟基的二醇单元 C 半缩醛+醚 SP1r 中等极性，环氧和羟甲基组合 D N-乙酰基 SP3d 酰胺官能团，极性 VS 虚拟位点 TC4 疏水珠子，无质量，增强π堆积相互作用葡萄糖醛酸（D-Glucuronic Acid，GlcA或GLA）化学结构：$\ce{C6H10O7}$（末端葡萄糖变为羧酸）映射原理：与葡萄糖类似，但C6($\ce{-CH2OH}$)被替换为羧基($\ce{-COOH}$) 原子级结构：C1-O1-C2($\ce{OH}$)-C3($\ce{OH}$)-C4($\ce{OH}$)-C5-O5-C6($\ce{COOH}$)，其中O1为异头体氧，O5为环氧（ether oxygen）粗粒化映射（四个珠子+虚拟位点）：珠子包含原子说明 A珠 C1-O1-C2 包含异头体碳C1和异头体氧O1 B珠 C3-C4 二醇单元 C珠 C5-O5 包含环氧O5 D珠 C6($\ce{COOH}$) 羧酸官能团，生理pH下去质子化 VS 虚拟位点 TC4类型，放置在环中心珠子类型选择依据：珠子碎片类型 Martini珠子类型选择依据 A 异头体 SP4r 异头体碳+O1，极性碎片 B 二醇 SP4r 标准二醇单元，含两个羟基 C 环氧醚 TN4ar 环氧和邻近碳 D 羧酸根 SQ5n（带电-1）生理pH下去质子化，强极性 VS 虚拟位点 TC4 增强π堆积相互作用实验分配系数验证（Table S2）：单糖实验Log P Martini 3预测（kJ/mol）误差（kJ/mol）精度评价 NAG -3.03 ± 0.34 -16.02 ± 0.33 1.27 优秀 GLA -3.26 ± 0.11 -18.17 ± 0.31 0.44 最优两种单糖的辛醇-水分配系数预测均达到高精度，验证了珠子类型选择和非键参数的准确性。内部环约束的15%缩放见正文Figure 1c，2（附录）。为了准确再现碳水化合物的分子体积和溶剂可及表面积（SASA），Martini 3对单糖环内的所有键长进行了统一的15%放大处理：环内键长：A-B、A-C、B-C（形成糖环的三个珠子之间的键）统一放大15% 糖苷键：连接两个单糖单元的键（如NAG的A珠到GlcA的B珠）不缩放，保持原始距离物理意义：直接从几何中心（COG）映射会低估分子体积约8%，15%的键长放大可使CG模型的Connolly表面与全原子参考高度一致适用性：这个缩放因子对所有单糖都适用，保证了模型的可迁移性单糖内部键合键合类型：使用约束（constraints）而非简谐键，因为单糖在CG层级表现为刚性三角形无angles/dihedrals：单糖环内三个珠子（A-B-C）之间不需要角度或二面角参数原始力场文件 [ moleculetype ] ; molname nrexcl GLA 1 [ atoms ] ; nr type resnr residue atom cgnr charge mass 1 SP4r 1 GLA A 1 0 54 2 SP4r 1 GLA B 2 0 54 3 TN4ar 1 GLA C 3 0 36 ; 4 SP3 1 GLA D 4 0 54 4 SQ5n 1 GLA D 4 -1.0 54 ;deprotonated at physiological pH 5 TC4 1 GLA VS 5 0 0 [constraints] ; i j funct length 1 2 1 0.376 ;15% COG scaled 1 3 1 0.335 2 3 1 0.311 3 4 1 0.222 ;unscaled, constraint because Fk > 80000 [angles] ; i j k funct angle fk 1 3 4 10 180 290 [dihedrals] ; i j k l funct angle fc 4 1 2 3 2 55 140 [ exclusions ] 5 1 2 3 4 4 2 [ virtual_sitesn ] 5 1 1 2 3 [ moleculetype ] ; molname nrexcl GYN 1 [ atoms ] ; nr type resnr residue atom cgnr charge mass 1 SN6 1 GYN A 1 0 54 2 SP4r 1 GYN B 2 0 54 3 SP1r 1 GYN C 3 0 54 4 SP3d 1 GYN D 4 0 54 5 TC4 1 GYN VS 5 0 0 [bonds] ; i j funct length fk 1 4 1 0.339 4700 ;unscaled [constraints] ; i j funct length 1 2 1 0.392 ;15% COG scaled 1 3 1 0.427 2 3 1 0.397 [ angles ] ; i j k funct angle fk 3 1 4 10 147 100 [dihedrals] ; i j k l funct angle fc 4 3 2 1 2 0 160 [ exclusions ] 5 1 2 3 4 4 2 [ virtual_sitesn ] 5 1 1 2 3 多糖建模图4：寡糖和多糖的参数化策略（详细讲解见下） a) 复杂碳水化合物的系统化映射策略 b) 两个连接的单糖片段之间引入的角度和二面角 c) 三个连续单糖片段之间引入的二面角 d) 糖苷键形成时新产生片段的珠子分配第一组（1-1、1-2、1-3、1-4链接）：使用SP1r珠子这个珠子类型直接来自单糖中的半缩醛片段已通过海藻糖和蔗糖的转移自由能验证（误差<3 kJ/mol）第二组（1-5、1-6链接）：使用SN6r珠子与半缩醛片段类似，但一个OH被醚键取代 SN6r的自相互作用比SP1r弱一级，反映了化学结构变化特殊情况（N-乙酰神经氨酸的1-4链接）：将羧酸与剩余碳片段组合，避免产生键长过短的2:1映射片段使用标准羧基珠子类型糖苷键参数化透明质酸（HA）的组成：由NAG（GlcNAc）和GlcA通过β-1,3糖苷键交替连接而成。糖苷键的分类 Martini 3将糖苷键分为六组，根据α/β异构体和链接碳位置：糖苷键类型例子映射方向接收方珠子类型 Class 1 α/β-1,1 & 1,2 异头体相连 T珠子 Class 2 α/β-1,3 & 1,4 最常见的β-1,4 T珠子 Class 3 α/β-1,5 & 1,6 包括6-脱氧 SN6r珠子（减弱相互作用）透明质酸中的β-1,3链接属于Class 2：这是该力场中最常见的链接类型之一。如何确定“接收单糖单元”？在糖苷键连接中，需要明确哪个单糖是“供体”（donor），哪个是“接收者”（acceptor）：规则：采用CHARMM-GUI约定，连接原子归属于CG层级中珠子编号更高的单糖单元例子：乳糖（α-1,4连接的葡萄糖-半乳糖）原子级连接：葡萄糖的C1连接到半乳糖的C4 CG级连接：葡萄糖的A珠连接到半乳糖的B珠糖苷醚氧原子归属于B珠（即半乳糖一侧，珠子编号更高的单元） β-1,3糖苷键的具体连接方式对于透明质酸的NAG-GlcA重复单元：原子级：NAG的C1（异头体碳）连接到GlcA的C3 CG级：NAG的A珠连接到GlcA的B珠糖苷醚氧归属：包含在GlcA的B珠中（接收方单糖）体积损失补偿糖苷缩合反应使总重原子数减少1（损失一个氧原子）：$\ce{C6H12O6 + C6H10O7 - H2O -> C12H20O11}$ Martini 3的解决方案：供体单糖（提供异头体碳C1的一侧）：保持原有珠子类型接收单糖（通过其他碳如C3/C4接收连接的一侧）：将接收糖苷键的珠子从S珠改为T珠（更小），以补偿重原子损失具体到透明质酸： NAG单元（供体）：A(SP1r) - B(SP1r) - C(SP1r) GlcA单元（接收方）：A’(TP1) - B’(SP1r，包含糖苷醚氧) - C’(SQ4) 注意：GlcA的A’珠从SP1r改为TP1（T珠），补偿糖苷缩合的重原子损失键合相互作用多糖键合参数糖苷键键长：从全原子参考映射获得，α和β异构体的键长明显不同，需分开处理 Angles（键角）：定义所有跨越两个单糖单元之间糖苷键的角度例如：A-糖苷键-B’，B-糖苷键-A’，A-糖苷键-C’等具体数值需从全原子MD模拟的分布拟合调和势获得 Dihedrals（二面角）：单糖内部：使用improper dihedral（funct=2，调和势）维持环平面性例如：GLA的4-1-2-3，用于保持糖环的平面构象主二面角（两个单糖连接）：使用proper dihedral（funct=1，周期性势函数）控制绕糖苷键的旋转（见Figure 4b）对于每个糖苷键，定义一个主二面角来控制绕该键的旋转二面角的具体原子选择取决于糖苷键连接类型（不同连接方式有不同的原子组合）例如：LAC (β-1,4链接，糖苷键为B-A’): 主二面角为A-B-A’-B’ 例如：SUCR/TREH (α-1,1链接，糖苷键为A-A’): 主二面角为B-A-A’-C’ 长程二面角（三个或更多单糖连接）：当连接超过两个单糖单元时，引入跨越三个连续单糖单元（n, n+1, n+2）的长程二面角，定义n和n+2残基相对于n+1残基平面的取向（见Figure 4c）对于含有N个单糖的多糖链，需要定义N-2个这样的长程二面角（每个连续三联体一个）例如：透明质酸（HA）的NAG₁-GlcA₂-NAG₃片段，长程二面角为B₁-A₂-B₂-A₃（从第1个残基选B珠，从第2个残基选A和B珠定义平面，从第3个残基选A珠），B₂-A₃-B₃-A₄，…… 这类二面角对多糖刚度至关重要，尤其是在较长的碳水化合物链中所有二面角参数通过匹配全原子参考模拟的构象分布获得受限弯曲势：对于被二面角势覆盖的角度，使用Bulacu等人的受限弯曲势，防止角度变为共线导致数值不稳定特殊处理葡聚糖（dextran）使用3-bonded neighbor exclusions以改善稳定性其他模型仅排除1-bonded neighbors（Martini脂质标准）虚拟位点的包含显著影响聚集行为和化学性质建模流程总览 mindmap root(碳水化合物建模) **单糖建模** 映射策略 **从C1逆时针分组** 最大化二醇单元保持官能团完整珠子分配基本3珠子：A-B-C A珠：异头体碳+O1 B珠：二醇单元 C珠：环氧O5 侧链D珠：NAG/GLA N-乙酰基：SP3d 羧基：SQ5n带电荷 **虚拟位点TC4**：π堆积键合参数 Constraints：环内键 **15%键长缩放** Improper dihedral：平面性 **多糖建模** 糖苷键规则 **糖苷醚氧归属珠子编号更高单元** **接收方S珠改为T珠**：补偿重原子损失 α/β键长不同需分开处理糖苷键分类 Class 1：α/β-1,1 & 1,2 Class 2：α/β-1,3 & 1,4 Class 3：α/β-1,5 & 1,6 键合参数糖苷键：不缩放 Angles：跨糖苷键角度主dihedral：单个糖苷键旋转 **长程dihedral：N-2个**，跨3残基参数化流程 1.全原子MD模拟 2.映射到CG珠子 3.拟合分布获参数 **验证与应用** 验证指标 SASA：小于5%误差转移自由能：1.5 kJ/mol **渗透压：解决粘性效应** 应用案例葡聚糖溶液性质蛋白质-糖脂识别糖蛋白/LPS体系验证方法与应用案例 Martini 3碳水化合物力场经过验证，在多个物理化学性质和实际应用中表现优异。详细内容请参见：附录：验证方法与应用案例验证指标概览力场验证基于三个核心物理化学性质：溶剂可及表面积（SASA） 15%键长缩放后，偏差 <5%（Martini 2为~8%） Connolly表面与全原子高度一致辛醇-水转移自由能平均绝对误差：1.5 kJ/mol 达到小分子Martini参数的精度水平渗透压 0-1.5 molal浓度：与实验优异吻合解决了Martini 2的”粘性效应”问题应用案例概览葡聚糖溶液性质：准确预测黏度、回转半径、扩散系数蛋白质-糖脂识别：成功模拟LecA与GM1的特异性结合糖蛋白、LPS、糖脂筏等复杂体系

Molecular Dynamics · 2025-11-16

Polyply：图匹配算法驱动的聚合物模拟参数生成与结构构建

Polyply：图匹配算法驱动的聚合物模拟参数生成与结构构建本文信息标题: Polyply; a python suite for facilitating simulations of macromolecules and nanomaterials 作者: Fabian Grünewald, Riccardo Alessandri, Peter C. Kroon, Luca Monticelli, Paulo C. T. Souza, Siewert J. Marrink 发表时间: 2022年1月单位: University of Groningen (荷兰格罗宁根大学) University of Chicago (美国芝加哥大学) CNRS and University of Lyon (法国里昂大学) 引用格式: Grünewald, F., Alessandri, R., Kroon, P. C., Monticelli, L., Souza, P. C. T., & Marrink, S. J. (2022). Polyply; a python suite for facilitating simulations of macromolecules and nanomaterials. Nature Communications, 13(1), 68. https://doi.org/10.1038/s41467-021-27627-4 GitHub代码: https://github.com/marrink-lab/polyply_1.0 文档: https://polyply.readthedocs.io Polyply官方文档 GitHub代码库 Martini力场官网摘要分子动力学模拟在（纳米）材料理性设计和生物大分子研究中扮演着日益重要的角色。然而，为这些模拟生成输入文件和真实的初始坐标是一个主要瓶颈，特别是对于高通量筛选协议和复杂多组分体系。为解决这一瓶颈，本文提出了Polyply软件套件，它提供：1）一个多尺度图匹配算法，能够快速生成任意复杂聚合物拓扑的参数；2）一个通用的多尺度随机游走协议，能够高效地设置复杂体系，且独立于目标力场或模型分辨率。作者通过创建聚合物熔体、单链及环状单链DNA的真实坐标来评估该方法的质量和性能，并通过设置微相分离嵌段共聚物体系和脂质囊泡内液-液相分离体系展示了该方法的强大功能。核心结论 Polyply基于图转换算法，将残基图（residue graph）映射为高分辨率参数文件，支持任意复杂的聚合物拓扑结构采用多尺度随机游走生成初始坐标，先构建超粗粒化（super CG）模型，再反向转换到目标分辨率力场无关设计，同时支持全原子和粗粒化模型，极大提升了高通量筛选的可行性在聚合物熔体、DNA单链、嵌段共聚物、相分离体系等多个复杂案例中验证了方法的准确性和效率背景分子动力学（MD）模拟已成为补充实验研究的强大工具。近年来，研究趋势从单一聚合物熔体或混合物转向更复杂的多组分体系，包括纯合成材料和生物-合成杂化大分子。这些体系的应用范围广泛，从聚电解质复合凝聚体到下一代聚合物电池，再到抗菌聚合物和可生物降解聚合物。随着材料基因组计划的推进，基于MD的虚拟高通量筛选正成为研究热点。MD高通量筛选相比实验方法成本更低，且能提供实验难以获取的性质信息，使研究者能够更高效地探索组合空间并筛选候选材料。然而，这一前景的实现需要程序能够快速、可靠、一致地构建拓扑和模拟盒子。当前的主要挑战在于：现有工具主要针对蛋白质、脂质膜、DNA等生物分子，对合成聚合物和生物-合成杂化大分子的支持严重不足。虽然存在一些特定解决方案，但它们通常只支持单一力场，仅限于开发者实现的特定（主要是线性）聚合物，且网站实现方案依赖服务器负载并需要人工交互。更复杂体系（如微相分离聚合物、杂化纳米颗粒共混物）的坐标生成往往依赖多尺度自组装或定制脚本。关键科学问题本文旨在解决聚合物和生物大分子MD模拟中的五个核心挑战：参数与坐标生成的通用性：程序需要同时生成坐标和参数，且与分辨率和力场无关。准确的粗粒化模型通常基于全原子聚合物，因此支持两者是高通量模型开发的关键输入文件生成的易用性：需要一个简单的流程，基于体系组成生成输入文件，支持任意复杂的聚合物序列，包括不同分支度和统计分布参数与坐标的组合能力：程序需要能够组合不同分辨率的聚合物输入，例如在相同模拟中混合全原子和粗粒化模型边界条件与几何的灵活性：需要支持三维周期性边界条件、球形、柱状等多种几何形状高通量筛选的性能要求：坐标和参数文件生成必须足够快，以支持高通量协议创新点图转换算法：首次将聚合物参数化问题完全转化为图同构匹配问题，实现了对任意复杂拓扑结构的自动参数生成多尺度随机游走：创新性地采用“超粗粒化→目标分辨率”的反向构建策略，避免了传统方法依赖坐标片段库的局限力场无关框架：通过分离算法核心与力场参数库，实现了对Martini、GROMOS、CHARMM、OPLS等多种力场的统一支持自动化工作流：从残基序列到完整模拟输入的全流程自动化，大幅降低了使用门槛研究内容 Polyply软件架构 Polyply由两个核心模块组成： polyply gen_params：基于图匹配算法的参数文件生成器 polyply gen_coords：基于多尺度随机游走的坐标生成器两个模块共享统一的图表示基础架构，均基于NetworkX和vermouth Python库实现图相关计算。图1：参数文件生成工作流程以聚乙二醇（PEO）接枝甲基丙烯酸酯（MA）为例，展示了三步图转换过程：输入：残基图（residue graph）和力场库中的building blocks 步骤1：生成目标分辨率的断开残基图步骤2：在残基图层级匹配links 步骤3：将通用links匹配到具体残基，生成完整参数文件核心算法一：图匹配驱动的参数生成 Polyply将参数文件生成问题转化为图转换（graph transformation）问题。其核心思想是：将残基图映射为高分辨率的分子图，该图与目标分辨率无关。基本概念图表示：分子的连接性转化为图的边，原子特征（名称、残基名等）存储为节点属性 Block（构建块）：对应单个残基的所有相互作用和原子的图 Link（连接）：描述两个或多个残基连接时引入的相互作用（如键、角度）三步图转换算法步骤1：生成断开的残基图遍历输入残基图的所有残基，为每个残基从库中匹配对应的block，添加到空图中，形成目标分辨率的断开图。此时已包含目标分子的所有原子和残基内相互作用，但缺少跨残基的连接。步骤2：在残基层级查找所有links 为生成跨越多个残基的相互作用，需要在残基之间应用links。Polyply将其转化为残基图层级的子图同构问题：查找link在残基图上的所有可能匹配方式，受节点属性等约束限制。在残基图层级执行大幅降低了问题规模。步骤3：将通用links匹配到具体残基根据步骤2建立的link与残基的对应关系，程序将link中的原子与步骤1生成的断开图中的原子建立对应关系。匹配不仅基于原子名称和残基索引，还可扩展到其他原子特征，从而考虑残基图连接性未编码的信息（如手性、端基异构体）。当link被添加时，其边也被添加到断开图中，逐步将断开图转变为目标分辨率的连通图。算法优势通用性：适用于任意复杂的聚合物拓扑，包括分支、环状、统计共聚等结构可扩展性：通过匹配节点属性，可处理手性、端基异构等精细化学信息效率：在残基图层级解决子图同构问题，显著降低计算复杂度核心算法二：多尺度随机游走坐标生成 Polyply采用通用多尺度方法构建起始坐标：首先生成超粗粒化（super CG）分辨率表示，然后反向转换到目标分辨率。这一策略类似于CHARMM-GUI polymer builder，但有三个关键改进：动态参数推导：super CG模型参数基于目标力场动态推导，而非预定义自排除随机游走：采用随机游走而非全尺度动力学模拟自动反向转换：不依赖坐标片段库的自动化反变换图2：坐标生成的五步工作流程五步坐标生成算法步骤1：将所有分子映射为每残基一个珠子分析拓扑文件，检测所有分子类型。对每个分子，识别所有唯一残基并转换为blocks。创建通用的每残基一个珠子的super CG模型，以图形式存储。残基图的连接性从分子的键合图中提取。步骤2：为残基生成坐标每个block是单个残基的图，使用图嵌入（graph embedding）生成坐标。由于分子几何的特殊要求，采用两步图嵌入：首先使用Kamada-Kawai嵌入生成初始坐标随后基于残基内键合相互作用进行几何优化，使用L-BFGS优化器步骤3：推导通用CG模型参数自排除随机游走使用每残基一个珠子的近似CG模型，基于Lennard-Jones（LJ）势。关键参数推导： ε参数（LJ势阱深度）：固定为1 kJ/mol（因不执行动力学，吸引部分不重要） σ参数（决定堆积密度）：从残基模板坐标计算，反映残基体积。基于回转半径推导（将聚合物物理中的链回转半径概念移植到单个残基的分子几何）此外，算法还考虑了残基在全原子模型中的天然堆积密度，通过缩放因子调整不同力场间的差异。步骤4：通过随机游走生成super CG坐标对体系中每个分子执行随机游走。算法依次添加残基：第一个残基随机放置后续残基通过以下方式添加：在前一个残基周围随机采样方向根据键合相互作用确定距离检查与已放置残基的重叠（使用LJ势）若无冲突则接受，否则重新采样这一过程确保了生成的构象满足键合约束，同时避免了原子重叠。步骤5：反向映射到目标分辨率将super CG坐标反向映射到目标分辨率。关键步骤：每个残基的质心固定在super CG珠子位置残基内部坐标从步骤2的模板继承应用适当的旋转和平移，确保跨残基键合几何正确对生成的结构进行能量最小化，消除局部应力坐标生成的关键技术多尺度策略：先在粗粒度生成全局构象，再细化局部结构，极大提升了效率自排除机制：随机游走过程中实时检测并避免原子重叠，确保生成结构的物理合理性自动反向映射：基于几何约束的自动化反变换，无需人工设计坐标片段库验证案例案例1：聚合物熔体作者测试了聚丙烯（PP）、聚乙烯（PE）、聚苯乙烯（PS）和聚甲基丙烯酸甲酯（PMMA）四种聚合物熔体的密度预测。结果：所有体系在5-10 ns内达到平衡密度误差<2%，与实验值高度一致 PP熔体（最苛刻测试）的Flory特征比与实验数据完美吻合这验证了Polyply生成的初始结构具有良好的物理性质，能快速弛豫到平衡态。案例2：单链DNA和环状DNA 作者使用Martini 3力场生成了单链DNA（ssDNA）和环状单链DNA（cssDNA）的坐标。图3：DNA结构生成与验证 a-c：ssDNA序列、生成的初始结构和平衡后的结构 d：ssDNA的末端距离分布与Martini 3全原子模拟高度一致 e-f：cssDNA的初始和平衡结构，展示了环状拓扑的正确处理关键发现：生成的ssDNA结构经短时间平衡后，末端距离分布与基准全原子模拟结果一致 cssDNA的环状拓扑约束得到正确处理，无需手动调整案例3：微相分离嵌段共聚物作者构建了聚苯乙烯-聚甲基丙烯酸甲酯（PS-PMMA）二嵌段共聚物的微相分离结构。图4：嵌段共聚物微相分离展示了PS-PMMA嵌段共聚物自组装形成的层状（lamellar）微相分离结构。图中不同颜色代表PS和PMMA嵌段，清晰显示了周期性层状相结构。结果： Polyply能够直接生成预组装的微相分离结构避免了耗时的自组装模拟过程生成的结构稳定，与已知相图一致案例4：脂质囊泡内的液-液相分离作者构建了一个复杂体系：脂质囊泡内包裹的液-液相分离（LLPS）体系。图5：脂质囊泡内的液-液相分离体系 a：体系组成示意图（脂质囊泡+LLPS液滴） b：生成的完整结构，展示了囊泡内两相分离的液滴技术亮点：演示了Polyply处理多组分、多尺度、复杂几何体系的能力组合了脂质（Martini粗粒化）、聚合物（LLPS相）、溶剂等多种组分支持球形约束等非周期边界条件性能评估图6：性能基准测试 a：参数生成时间随聚合物长度的缩放关系（线性缩放） b：坐标生成时间随聚合物长度的缩放关系 c：坐标生成成功率随体积分数的变化关键结论：参数生成对数千个残基的聚合物仅需秒级时间坐标生成时间随链长近似线性增长在高体积分数（φ > 0.5）下仍能保持>90%的成功率 Q&A Q1：Polyply的图匹配算法与传统参数生成方法相比有何优势？ A1：传统方法通常针对特定聚合物类型编写专门代码，扩展性差。Polyply的图匹配算法将问题抽象为通用的子图同构匹配，只需定义building blocks和links即可支持新聚合物类型，无需修改核心代码。此外，在残基图层级执行匹配大幅降低了计算复杂度。 Q2：多尺度随机游走为什么不直接在目标分辨率生成坐标？ A2：直接在目标分辨率（特别是全原子）执行随机游走面临巨大的构象空间采样问题，且容易产生原子重叠。先在super CG层级生成全局构象可以：1）大幅减少自由度，提升采样效率；2）更容易满足键合约束；3）通过LJ势简单有效地避免大尺度重叠。反向映射步骤则利用局部几何优化解决精细尺度的冲突。 Q3：Polyply如何确保生成的聚合物链长分布符合实验？ A3：Polyply允许用户指定任意的链长分布（单分散、多分散、特定分子量分布等）。用户可以通过输入文件定义每条链的确切序列，或使用统计分布函数（如高斯分布、指数分布）来模拟真实的分子量分布。这为模拟真实聚合物样品提供了灵活性。 Q4：对于高度分支的聚合物（如树枝状大分子），Polyply是否适用？ A4：是的。Polyply的图表示天然支持任意拓扑结构，包括高度分支、星形、树枝状等。只需在残基图中正确定义分支点的连接关系，算法会自动处理所有跨残基的相互作用。作者在文中已演示了接枝共聚物（PEO-g-MA）的参数生成。 Q5：Polyply生成的初始结构质量如何？是否需要长时间平衡？ A5：从基准测试来看，Polyply生成的结构质量很高。聚合物熔体案例中，体系在5-10 ns内即达到平衡密度；DNA案例中，末端距离分布经短时间平衡后与全原子基准一致。这表明生成的结构已接近物理合理的构象，大大缩短了后续模拟的平衡时间。关键结论与批判性总结主要贡献 Polyply通过图转换算法实现了聚合物参数化的完全自动化，支持任意复杂拓扑结构多尺度随机游走策略在保证坐标质量的同时显著提升了生成效率力场无关的软件架构使其能广泛应用于不同力场和模型分辨率在聚合物熔体、DNA、嵌段共聚物、LLPS等多个复杂体系的成功应用验证了方法的鲁棒性局限性高体积分数限制：虽然在φ > 0.5时仍有>90%成功率，但对于极高密度体系（如晶体），随机游走方法可能需要过多尝试力场库依赖：虽然用户可自定义blocks和links，但对于全新化学体系，仍需手动构建参数库环状聚合物的闭环约束：对于大环聚合物，反向映射后闭环可能引入较大应力，需要更仔细的能量最小化动力学性质：论文主要验证了结构和热力学性质，对于依赖精确动力学的性质（如扩散系数、粘度）的适用性需进一步验证未来展望参数库扩展：建立涵盖更多化学单元的社区参数库，提升开箱即用性机器学习集成：利用ML预测最优super CG参数，进一步提升坐标生成效率晶体结构支持：开发针对晶格结构的专门算法，扩展到聚合物晶体模拟与实验数据整合：结合散射实验数据（SAXS、SANS）优化生成结构，提升与实验的一致性总体评价 Polyply代表了聚合物模拟工作流自动化的重大进步。其通用的图算法框架和力场无关设计，使其能够成为连接不同力场、不同分辨率、不同聚合物类型的统一平台。特别是对于高通量虚拟筛选这一新兴应用，Polyply提供的快速、自动化工作流具有不可替代的价值。虽然仍存在一些局限性，但软件的开源性和模块化设计为社区贡献和持续改进提供了良好基础。

Molecular Dynamics · 2025-11-16

Polyply技术细节：算法实现与扩展案例（附录）

本文是《Polyply：图匹配算法驱动的聚合物模拟参数生成与结构构建》的附录，包含详细的算法实现、参数推导和扩展验证案例。算法实现细节图嵌入与几何优化 Polyply使用两步图嵌入策略生成残基的初始坐标：步骤1：Kamada-Kawai嵌入 Kamada-Kawai算法将图嵌入问题转化为能量最小化： \[E = \sum_{i<j} k_{ij} (d_{ij} - l_{ij})^2\] 其中： $d_{ij}$是节点i和j之间的欧几里得距离 $l_{ij}$是图中i和j之间的最短路径长度 $k_{ij} = K / l_{ij}^2$是弹簧常数该算法能生成反映图拓扑的初始坐标，但不考虑分子几何约束。步骤2：L-BFGS几何优化基于残基内键合相互作用进行几何优化，目标函数： \[F = \sum_{\text{bonds}} k_b (r - r_0)^2 + \sum_{\text{angles}} k_\theta (\theta - \theta_0)^2 + \sum_{\text{dihedrals}} k_\phi [1 + \cos(n\phi - \delta)]\] 使用L-BFGS算法最小化，确保生成的残基几何满足力场约束。 Super CG模型参数推导回转半径计算对于单个残基，回转半径定义为： \[R_g = \sqrt{\frac{1}{N} \sum_{i=1}^{N} (\mathbf{r}_i - \mathbf{r}_{\text{COM}})^2}\] 其中$\mathbf{r}_{\text{COM}}$是质心坐标。 LJ σ参数推导 super CG模型的σ参数基于回转半径： \[\sigma = 2 R_g \times f_{\text{scale}}\] 缩放因子$f_{\text{scale}}$根据力场调整： GROMOS全原子：$f_{\text{scale}} = 1.0$ Martini粗粒化：$f_{\text{scale}} = 0.85$ 这一差异反映了不同力场中残基天然堆积密度的不同。自排除随机游走算法伪代码如下：对于每个分子：将第一个残基随机放置在盒子中对于后续每个残基： max_attempts = 1000 for attempt in range(max_attempts): # 随机采样方向 direction = random_unit_vector() # 根据键长确定距离 distance = bond_length(previous_residue, current_residue) # 计算候选位置 candidate_position = previous_position + distance * direction # 检查与所有已放置残基的重叠 overlap = False for placed_residue in placed_residues: LJ_energy = calculate_LJ(candidate_position, placed_residue) if LJ_energy > threshold: # 默认10 kJ/mol overlap = True break if not overlap: accept_position(candidate_position) break if overlap: # 所有尝试都失败 return FAILURE 关键参数：重叠阈值：10 kJ/mol（对应约0.7σ的距离）最大尝试次数：1000次/残基扩展验证案例聚合物熔体详细数据作者测试了多种聚合物熔体，详细数据见下表：聚合物力场温度(K) 实验密度(g/cm³) 模拟密度(g/cm³) 误差(%) PP GROMOS 513 0.76 0.74 ± 0.01 2.6 PE GROMOS 413 0.78 0.77 ± 0.01 1.3 PS GROMOS 513 0.97 0.95 ± 0.02 2.1 PMMA GROMOS 513 1.10 1.08 ± 0.02 1.8 PEO Martini 413 1.06 1.05 ± 0.01 0.9 PMA Martini 413 1.10 1.09 ± 0.01 0.9 所有体系在5-10 ns内达到平衡密度，表明Polyply生成的初始结构质量高。 DNA末端距离分布 SI图1：DNA回转半径和末端距离分布左图：回转半径分布右图：末端距离分布红色：全原子MD模拟参考蓝色：Polyply生成的200个初始结构关键观察： Polyply生成的分布较宽，但与全原子分布有良好重叠证明Polyply构象是良好的起始点注意：全原子力场预测的$R_g = 2.8 \pm 0.5$ nm低于实验值$3.8 \pm 0.1$ nm 环状DNA在病毒衣壳内的构建 SI图2：猪病毒环状ssDNA生成工作流程案例亮点：从数据库获取病毒基因组序列和衣壳晶体结构使用ParmSC1力场为DNA生成参数衣壳蛋白使用Amber14力场 DNA使用球形几何约束+衣壳边界限制每个核苷酸位点放置一个$\ce{Na+}$离子（使用ligation功能）使用cycle选项生成环状DNA 关键技术：球形约束加速算法（避免与每个衣壳原子检查重叠）高盐浓度（~2 mol/L）使DNA采用柔性无规卷曲构象三步平衡：0.1 fs柔性键 → 1 fs约束键 → 2 fs生产运行结果：60 ns生产运行中体系稳定，观察到衣壳内外的离子交换，暗示衣壳内存在最优盐浓度。聚合物锂离子电池 SI图3：PS-b-PEO LiTFSI掺杂电池生成工作流程体系组成：聚苯乙烯-聚乙二醇二嵌段共聚物（PS-b-PEO）锂双三氟甲烷磺酰亚胺盐（LiTFSI）掺杂 Martini 2粗粒化力场验证结果：层间距：模拟值~21 nm，实验值20 nm（优异吻合）盐分布：$\ce{Li+}$富集在PEO畴内，与实验报道的盐通道形成一致相分离：清晰的PS和PEO交替层状结构，界面有一定混合这一案例展示了Polyply在功能材料模拟中的应用潜力。脂质囊泡内液-液相分离详细工作流程 SI图4：葡聚糖-PEO液-液相分离囊泡工作流程葡聚糖分子量分布建模：作者使用线性缩聚反应动力学模型： \[\text{prob}(N, p) = N \times p^{N-1} (1-p)^2\] 其中$p$是反应程度。通过调整$p$使数均分子量$\bar{M}_n \approx 65$（与实验一致），得到多分散指数PDI $\approx 1.5$（文献值1.8）。支化度：5%的1,3-糖苷键（文献值，分子量<100,000 g/mol）结果：成功生成包含500个不同链长葡聚糖分子的多分散体系展示了Polyply处理统计共聚和多分散性的能力性能优化策略参数生成优化子图同构匹配：在残基图层级执行而非原子层级，复杂度从$O(N_{\text{atoms}}!)$降至$O(N_{\text{residues}}!)$ 缓存机制：相同残基类型的block只需加载一次并行化：独立分子的参数生成可并行执行坐标生成优化 Early termination：检测到不可能完成的构象立即终止（如体积分数过高）分层放置：优先放置大分子，小分子填充空隙网格加速：使用空间分区网格加速重叠检测，复杂度从$O(N^2)$降至$O(N \log N)$ 成功率与体积分数作者系统测试了不同体积分数下的成功率：体积分数φ 成功率平均尝试次数/残基 0.1 100% <10 0.3 99% <50 0.5 95% <200 0.7 90% <500 0.9 <50% >1000 建议： φ < 0.7：直接使用Polyply 0.7 < φ < 0.9：增加max_attempts或使用更小的初始盒子尺寸 φ > 0.9：考虑先在较低密度生成，再通过NPT压缩力场库扩展当前支持的力场全原子：GROMOS 54A7, GROMOS 2016H66, Amber14, CHARMM36 粗粒化：Martini 2, Martini 3, SDK（软球模型）添加新残基示例创建一个PEO单元的block文件（JSON格式）： { "name": "PEO", "atoms": [ {"name": "C1", "type": "CH2", "charge": 0.0}, {"name": "O", "type": "O", "charge": -0.4}, {"name": "C2", "type": "CH2", "charge": 0.0} ], "bonds": [ {"atoms": ["C1", "O"], "length": 0.143, "force_constant": 8000}, {"atoms": ["O", "C2"], "length": 0.143, "force_constant": 8000} ], "angles": [ {"atoms": ["C1", "O", "C2"], "angle": 109.5, "force_constant": 450} ] } 创建对应的link文件定义C2-C1’连接： { "name": "PEO-PEO", "atoms": ["C2", "+C1"], "bond": {"length": 0.153, "force_constant": 7500} } 常见问题与解决方案问题1：坐标生成失败症状：生成过程卡住或报错“Maximum attempts reached” 可能原因：体积分数过高残基间存在不兼容的几何约束 LJ参数设置不合理解决方案：降低目标密度，稍后通过NPT压缩检查残基模板坐标的合理性调整$f_{\text{scale}}$参数问题2：生成结构需要长时间平衡症状：能量最小化或MD平衡耗时过长可能原因：存在严重的原子重叠键长/键角与力场参数偏差大解决方案：降低重叠阈值（更严格的重叠检测）使用更精细的几何优化（增加优化步数）分阶段平衡（逐步增加时间步长）问题3：环状聚合物闭环失败症状：环不闭合或闭环处应力过大可能原因：链长与持久长度不匹配随机游走未考虑闭环约束解决方案：使用更灵活的链（降低持久长度）先生成开链，后通过约束MD逐步闭合增加Monte Carlo尝试次数与其他工具的比较特性 Polyply CHARMM-GUI Packmol Moltemplate 参数生成 ✓ ✓ ✗ ✓ 坐标生成 ✓ ✓ ✓ ✗ 力场无关 ✓ ✗ ✓ ✓ 任意拓扑 ✓ 部分 ✗ ✓ 高通量友好 ✓ ✗ ✓ 部分图形界面 ✗ ✓ ✗ ✗ Polyply的独特优势：唯一同时支持参数和坐标生成、且力场无关的工具图算法框架提供最大的灵活性和可扩展性命令行界面最适合高通量脚本化工作流未来技术路线图机器学习增强：使用ML预测最优super CG参数和重叠阈值 GPU加速：将重叠检测和能量计算移至GPU 云服务：提供Web界面和REST API，降低使用门槛与自动化力场开发工具集成：如GAFF、CGenFF自动参数化工具晶格结构模板：为聚合物晶体提供专门的构建算法相关资源主文档：Polyply核心原理和主要应用 Polyply官方教程 GitHub Issues：问题反馈和讨论

Molecular Dynamics · 2025-11-16

Martini 3 脂质组学：更精细的参数如何重塑膜模拟的未来

Martini 3 脂质组学：更精细的参数如何重塑膜模拟的未来本文信息标题: Martini 3 脂质组学：扩展和精炼的参数改善脂质相行为作者: Kasper B. Pedersen, Helgi I. Ingólfsson, Siewert J. Marrink, Paulo C. T. Souza 等 (多国合作团队) 发表时间: 2025年7月31日单位: 奥胡斯大学 (丹麦)，劳伦斯利弗莫尔国家实验室 (美国)，卡尔加里大学 (加拿大)，格罗宁根大学 (荷兰) 等引用格式: Pedersen, K. B., Ingólfsson, H. I., Ramirez-Echemendia, D. P., Borges-Araújo, L., Andreasen, M. D., Empereur-mot, C., … & Marrink, S. J. (2025). The Martini 3 Lipidome: Expanded and Refined Parameters Improve Lipid Phase Behavior. ACS Central Science, 11, 1598–1610. https://doi.org/10.1021/acscentsci.5c00755 源代码/数据库: https://github.com/Martini-Force-Field-Initiative/M3-Lipid-Parameters 摘要脂质膜是细胞生命的核心。作为实验的补充，计算模拟在揭示复杂的脂质-生物分子相互作用方面至关重要，无论在学术界还是工业界都扮演着关键角色。Martini模型，一种用于高效分子动力学模拟的粗粒化力场，被广泛用于研究膜现象，但也面临着局限性，特别是在捕捉真实的脂质相行为方面。在这里，我们提出了一套精炼的Martini 3脂质模型，其采用的映射方案能够区分仅相差两个碳原子的脂质尾链，从而增强了包括三元混合物在内的模型膜系统的结构分辨率和热力学准确性。扩展后的Martini脂质库包含了数千个模型，使得对复杂且具有生物学相关性的系统进行模拟成为可能。这些进展将Martini确立为一个跨越多个领域的、强大的脂质模拟平台。核心结论提出了全新的Martini 3脂质映射方案：通过引入小尺寸珠子，新方案能够区分长度仅相差2个碳原子的脂质尾链（例如16C vs 18C），极大地提升了模型的化学分辨率。构建了庞大的脂质库：通过自动化脚本和精细的参数化流程，生成了包含数千种不同脂质的Martini 3模型库，涵盖了多种头基和尾链组合。显著改善了相行为的预测：与Martini 2相比，新的Martini 3脂质模型在预测脂质的凝胶-液晶相变温度 ($T_m$) 和三元混合物（如DPPC/DOPC/CHOL）的液有序(Lo)/液无序(Ld)相分离方面，与实验数据达到了前所未有的吻合度。提升了膜力学性质的准确性：新模型计算出的膜弯曲模量 ($k_c$) 和脂质尾链有序度参数也比Martini 2更接近全原子模拟的结果。成功模拟了复杂生物膜与非层状结构：展示了新脂质组学在构建真实的、不对称的哺乳动物细胞质膜模型以及模拟反相六方相和立方相等对药物递送至关重要的非层状结构中的强大能力。背景细胞膜是生命活动的基础舞台，它不仅是细胞的物理边界，更是无数生物化学反应发生的场所。从蛋白质折叠到信号转导，再到病毒入侵，几乎所有关键生命过程都与膜的结构和动态特性息息相关。然而，膜的复杂性——由成百上千种不同的脂质分子动态组成——使得单纯的实验研究难以捕捉其全貌。因此，分子动力学 (MD) 模拟，特别是粗粒化 (Coarse-Grained, CG) 模拟，已成为膜生物物理学研究不可或缺的工具。在众多CG模型中，Martini力场以其高效与准确的平衡而独树一帜，成为过去二十年中最流行的CG力场之一。它通过将多个原子“打包”成一个相互作用珠子，极大地降低了计算复杂度，使得模拟的时间和空间尺度可以达到微秒和数百纳米级别，从而能够研究膜的自组装、相分离（脂筏的形成）和与蛋白质的相互作用等宏观现象。然而，尽管Martini 2版本取得了巨大成功，但它也存在着一些众所周知的局限。其中最突出的一个便是对脂质相行为的描述不够准确。例如，Martini 2的映射方案无法区分DPPC (16:0) 和DSPC (18:0)这两种饱和脂质，尽管它们的相变温度在实验中相差14度之多。更重要的是，在模拟经典的DPPC/DOPC/胆固醇三元混合物时，Martini 2无法重现实验中观察到的液有序(Lo)-液无序(Ld)相分离，这极大地限制了其在研究细胞膜上功能性微区（如脂筏）时的可靠性。随着Martini 3的发布，其更丰富的珠子类型和更灵活的参数化策略为解决这些难题提供了契机。关键科学问题本文旨在对Martini 3的脂质模型进行一次系统性、大规模的重参数化和扩展，以解决Martini 2的上述局限性。其核心科学问题可以分解为：如何提高模型的化学分辨率？能否设计一种新的映射方案，使其能够精确地区分化学结构上仅有细微差异（如相差两个亚甲基）的脂质分子？如何平衡“自下而上”与“自上而下”的参数化策略？能否开发一套流程，既能保证CG模型在局部结构（如键长、键角分布）上与全原子模拟（“自下而上”）匹配，又能确保其宏观性质（如膜厚、相变温度）与实验数据（“自上而下”）吻合？新模型是否真正解决了核心痛点？经过重新参数化后，新的Martini 3脂质模型在预测凝胶-液晶相变和三元体系相分离这两个经典难题上的表现究竟如何？新模型的适用性有多广？这套经过优化的参数和模型是否能够被推广，用于构建包含数十种脂质的真实生物膜模型，并准确模拟其与蛋白质的相互作用以及非层状相的形成？创新点全新的脂质映射方案：创造性地引入了“小尺寸”珠子，实现了对脂质尾链长度每2个碳原子进行区分的能力，彻底解决了Martini 2中因“模糊映射”导致的不同脂质共用同一模型的问题。混合式参数化工作流：建立了一套严谨的“两阶段”参数化流程。第一阶段，通过与CHARMM36全原子模拟的键长、键角分布进行拟合，确保局部结构的准确性；第二阶段，通过与一个大型实验数据库（本文称之为MIB）中的宏观性质（如膜厚、相变温度）进行比对，进行“人工在环”的微调，确保了全局性质的真实性。建立了“Martini脂质基准” (MIB)：通过广泛的文献调研，整理并建立了一个包含29种脂质、在不同温度下共计67个数据点的公开实验数据库，为当前和未来的力场开发提供了一个宝贵的“黄金标准”。实现了数千种脂质的自动化建模：开发了一套自动化脚本，可以根据新的参数化构建块，快速生成数千种不同磷脂、鞘磷脂、神经酰胺等脂质的Martini 3拓扑文件，极大地扩展了Martini脂质组学。研究内容核心方法：两阶段参数化与实验基准验证本文的核心方法是一套结合了“自下而上”的精确性和“自上而下”的真实性的混合参数化策略。 graph TD subgraph "方向：从左到右" direction LR A["1.定义新的映射方案 引入小尺寸珠子 区分2个碳原子差异"] --> B["2.自下而上参数化 (Bottom-up) 构建全原子参考体系(CHARMM36) 拟合CG模型的键长、键角分布"]; B --> C["3.自上而下验证 (Top-down) 构建大型实验数据库(MIB) 模拟大量单组分膜体系"]; C --> D{"4.比较模拟与实验 (膜厚、APL、相变温度等)"}; D -- "不匹配" --> E["5.人工在环优化 (Human-in-the-loop) 微调参数以权衡各项性质"]; E --> B; D -- "匹配" --> F["最终优化的 Martini 3脂质参数"]; end 图1：重现各种主要脂质类别的结构性双层膜性质。 (A) Martini 3脂质模型的重新定义映射方案。(B-C) 参数化策略首先匹配高分辨率CHARMM36脂质模型的键和角分布，然后测试一系列涌现的双层膜性质，如双层膜几何形状和相行为。(D-G) 将Martini 3双层膜的模拟结果与大型实验基准（MIB）进行比较。 1. 全新的映射方案：更高的化学分辨率 Martini 2最大的问题之一是其“4对1”的映射规则过于粗糙。为了解决这个问题，作者在Martini 3的框架下引入了小尺寸珠子(S)。例如，对于一个16碳的棕榈酸链，Martini 3现在将其映射为 SNda-SC1-C1-C1 (一个S珠子加三个C珠子)，而对于18碳的硬脂酸链，则映射为 SNda-C1-C1-C1-C1 (四个C珠子)。这种精细的划分是实现对不同脂质精确描述的基础。 2. 两阶段参数化：从原子到宏观阶段一：自下而上拟合：研究人员首先进行了大量不同脂质的全原子模拟（使用CHARMM36力场），然后将这些轨迹“映射”成粗粒化的伪轨迹。接着，他们调整Martini 3的键长和键角参数，使得CG模拟的键长、键角分布函数与全原子伪轨迹的分布函数尽可能吻合。阶段二：自上而下验证：这是最关键的一步。作者整理了一个包含29种脂质在不同温度下的面积、厚度等实验数据的大型基准数据库（MIB）。他们用第一阶段得到的参数进行大量单组分膜的CG模拟，计算相应的宏观性质，并与MIB中的实验值进行比较。如果存在偏差，他们会“人工在环”地微调一些关键参数（例如饱和脂质尾链的角度力常数），在局部结构准确性和宏观性质真实性之间寻找最佳平衡点。结果与分析 1. 宏观性质与实验高度吻合图1 (D-G) 展示了新参数化后的Martini 3模型在预测四个关键宏观性质上的表现：单脂质面积 (APL): 模拟值与实验值的相关性极高，尽管为了改善相变行为，PC、PG和SM脂质的APL被有意地略微低估了约3 Å²。膜厚 (DHH, DB, 2Dc): 无论是磷酸头基间的峰-峰距离(DHH)、总厚度(DB)还是疏水核心厚度(2Dc)，模拟值都与实验值表现出极好的一致性。这些结果证明，新的映射方案和参数化策略成功地捕捉了不同脂质在形成双层膜时的几何特征。 2. 力学与动态性质的改进图S1：重新参数化的Martini 3脂质改善了双层膜弯曲模量和脂质尾链有序度。 (A) M2（旧版）和M3（新版）计算的弯曲模量($k_c$)与CHARMM36结果的比较。(B) M2和M3的$P_2$有序度参数与CHARMM36结果的比较。弯曲模量($k_c$)：这是一个描述膜抵抗弯曲能力的力学性质。如图S1A所示，新Martini 3模型（M3）计算的$k_c$值与全原子模拟（C36）的相关性（$R^2=0.97$）远高于旧的Martini 2模型（M2）。尾链有序度($P_2$)：如图S1B所示，对于多种脂质，M3的尾链有序度曲线（红色）也比M2（蓝色）更贴近全原子模拟（黑色）的结果。 3. 核心突破：精确预测相变温度这是本文最核心的成果之一。作者使用了两种互补的方法来确定相变温度 ($T_m$)。图2：改进的Martini 3脂质的相行为。 (A,C) 模拟退火方法。(B,D) 晶种法。模拟退火：通过缓慢地降低和升高温度，观察体系的有序度指标（林德曼指数）发生突变的位置（图2C）。晶种法 (Seeding)：在一个模拟盒子中同时放入一块固相（凝胶相）和一块液相的膜，然后在一系列不同温度下进行模拟，观察哪个相会“吞噬”另一个相，从而精确地“夹逼”出相变温度（图2D）。表1：饱和脂质的相变温度（开尔文）脂质模拟退火 $T_m$ 晶种法 $T_m$ 实验 $T_m$ DPPC (16:0/16:0) 323.1±3.5 320.0±4.0 314 DSPC (18:0/18:0) 338.5±3.5 330.5±2.5 328 PSM (d18:1/16:0) 323.8±3.0 313.5±2.5 314 SSM (d18:1/18:0) 322.5±1.5 314.0±1.0 318 结果令人振奋：新Martini 3模型预测的$T_m$值与实验值的误差在5 K以内，并且能准确地区分DPPC和DSPC。 4. 核心突破：重现三元体系相分离这是本文最核心的突破之一。作者采用了大规模三元混合物模拟的方法来判断DPPC/DOPC/胆固醇体系的相分离行为。模拟设置与技术细节：构建DPPC/DOPC/CHOL三元体系，在40 × 40 × 10 nm的模拟盒子中包含约6000个脂质分子和约130,000个总粒子。在297 K温度下进行10 μs的超长时间模拟，确保体系达到平衡并观察到相分离现象。关键技术参数：使用z轴位置约束：对上层膜的PO4珠子施加2 kJ/mol/nm²的位置约束，防止大尺度膜起伏影响相分离行为采用半各向同性压力耦合：在x-y平面内允许膜自由调整尺寸，同时保持z方向独立控制设置合适的邻居列表截断距离：1.35 nm，确保正确的邻居列表更新相分离判断与定量分析方法：视觉识别：通过分子动力学轨迹的快照，直接观察不同脂质组分是否形成宏观分离的区域。Lo相（富含DPPC和CHOL）和Ld相（富含DOPC）会在膜平面上形成清晰的相分离图案脂质富集分析：使用LiPyphilic等分析工具的Neighbours模块计算脂质邻居富集指数（Enrichment Index），定量描述不同脂质之间的聚集程度定量表征指标：密度分布曲线：计算不同组分沿膜法线方向的密度分布，Lo相和Ld相具有不同的脂质头基和尾链分布特征膜厚差异：Lo相由于DPPC和CHOL的紧密堆积，通常比Ld相具有更大的膜厚有序度参数：通过计算脂质尾链的P2有序度参数，Lo相显示更高的有序度值关键指标详解林德曼指数（Lindemann Index）林德曼指数是用于判断脂质尾链是否处于凝胶相的关键定量指标。它源于固体物理学，用来描述原子或分子在其平衡位置附近的均方根波动。计算公式：对于每个脂质尾链珠子i，在时间窗口内的林德曼指数定义为： \[\delta_i = \frac{1}{N-1} \sum_{j \neq i} \frac{\sqrt{\langle r_{ij}^2 \rangle - \langle r_{ij} \rangle^2}}{\langle r_{ij} \rangle}\] 其中： $r_{ij}$ 是珠子i与相邻珠子j之间的距离 $N$ 是尾链中的珠子总数 $\langle \cdots \rangle$ 表示时间平均该指数对所有相邻珠子对进行平均物理意义： $L < 0.1$：脂质尾链高度有序，处于凝胶相（Lβ）或固相，分子排列紧密，热运动受限 $0.1 < L < 0.15$：过渡区域，可能是液有序相（Lo）或接近相变点 $L > 0.15$：脂质尾链无序，处于液晶相（Ld），分子运动自由在本研究中的应用：在模拟退火实验（图2A和2C）中，作者通过监测林德曼指数随温度的变化来识别相变温度 $T_m$ 当体系从高温降温时，林德曼指数会在相变点附近发生突变（从>0.15突降至<0.1）这种突变对应于从液晶相到凝胶相的转变，其转折点即为相变温度富集指数（Enrichment Index）富集指数是用于定量描述三元混合物中脂质相分离程度的核心指标。它衡量某种脂质分子周围出现另一种脂质分子的概率是否偏离随机分布。计算方法：对于脂质类型A和B，富集指数 $E_{AB}$ 定义为： \[E_{AB} = \frac{N_{AB}^{\text{obs}}}{N_{AB}^{\text{exp}}} - 1\] 其中： $N_{AB}^{\text{obs}}$ 是实际观察到的A分子周围B分子的数量（通常定义为第一壳层内，如4-6 Å范围） $N_{AB}^{\text{exp}}$ 是基于随机分布预期的B分子数量，计算为：$N_{AB}^{\text{exp}} = N_{\text{total}} \times \frac{n_B}{n_A + n_B}$ 物理意义： $E_{AB} > 0$：A和B倾向于聚集在一起，表明两者相互吸引或倾向于共存于同一相 $E_{AB} = 0$：A和B的分布是随机的，不存在相分离 $E_{AB} < 0$：A和B倾向于分离，表明两者相互排斥或存在相分离在本研究中的应用：在图2H中，作者展示了Ca²⁺存在下POPS:POPC 50:50体系的富集指数图中显示POPC（蓝色线）和POPS（红色线）的富集指数： POPC-POPC富集指数显著为正：说明POPC分子倾向于聚集在一起，形成富POPC的Ld相 POPS-POPS富集指数显著为正：说明POPS分子（尤其是在Ca²⁺作用下）也倾向于聚集，形成富POPS的凝胶相区域 POPC-POPS交叉富集指数为负：说明两种脂质倾向于分离，证实了相分离的存在在DPPC/DOPC/CHOL三元体系的研究中，通过计算不同组分之间的富集指数，可以定量确认Lo相和Ld相的形成及其边界图2：改进的Martini 3脂质的相行为。 (A,C) 模拟退火方法。(B,D) 晶种法。(E-H) 离子诱导的相变：(E) Ca²⁺存在下POPS双层膜的液相→凝胶相转变；(F) Ca²⁺存在下POPS:POPC 50:50双层膜的相分离；(G) 从POPS双层膜模拟中获得的离子（Na⁺, Ca²⁺）、磷酸盐（PO₄）和甘油连接基（GLs）的密度分布；(H) 从Ca²⁺存在下POPS:POPC 50:50双层膜模拟中获得的POPC（蓝色）和POPS（红色）脂质的富集指数。(I-K) DPPC、DOPC和CHOL混合物的三元相行为：(I) 来自实验的相图；(J) 使用Martini 2模拟的相图；(K) 使用Martini 3模拟的相图。模拟结果与实验对比：实验相图 (图2I)：相图非常复杂，存在单相区（Ld, Lo）和多相共存区（Ld/Lo, Ld/Lβ, Lo/Lβ, Ld/Lo/Lβ） Martini 2相图 (图2J)：模拟结果非常糟糕，几乎整个相图都是均一的液无序相(Ld)，完全无法捕捉到相分离 Martini 3相图 (图2K)：模拟结果与实验惊人地吻合。不仅纯DPPC形成了正确的凝胶相(Lβ)，而且在正确的组分区域出现了Ld/Lβ和Lo/Lβ的相分离，甚至还捕捉到了一部分三相共存的区域三元相图解读三元相图（图2I-K）采用了蜂窝状六边形网格，每个六边形代表一个特定的DPPC/DOPC/CHOL组分比例，通过不同的颜色编码来表示该组分下的相态：单相区域：红色：纯液无序相（Ld），主要出现在高DOPC含量区域。特征是脂质尾链无序、膜较薄、流动性高绿色：纯液有序相（Lo），主要出现在高DPPC和高CHOL区域。特征是脂质尾链有序、膜较厚、但仍保持侧向流动性深紫色/黑色：纯凝胶相（Lβ），主要出现在高DPPC、低CHOL区域（CHOL浓度<20%）。特征是脂质尾链高度有序、膜最厚、侧向扩散几乎冻结两相共存区域：黄色/橙色：Ld + Lo相共存，这是最重要的生物学相关区域，对应于细胞膜上的”脂筏”现象。膜表面同时存在流动的无序区（富DOPC）和有序的微区（富DPPC+CHOL）蓝色/青色：Lo + Lβ相共存，常见于低CHOL、中等DPPC含量区域。膜表面同时存在流动相和凝胶相的岛屿粉色/浅紫色：Ld + Lβ相共存，出现在高DPPC、中等CHOL含量区域三相共存区域：白色或灰色：Ld + Lo + Lβ三相共存，这是相图中最复杂的区域，三种相态同时存在。只在非常窄的组分范围内出现关键发现对比：从图2的三个相图（I实验、J-M2、K-M3）对比可以看出：实验相图（I）的主要特征：左下角（高DOPC）为红色Ld相右下角（高DPPC，低CHOL）为粉色/浅紫色Lβ相右上角（高DPPC+高CHOL）为绿色Lo相存在明显的黄色Ld/Lo共存带、蓝色Lo/Lβ共存带和深绿色Ld/Lβ共存带 Martini 2的失败（J）：几乎整个相图都是红色（Ld相），只有最右下角极小区域显示凝胶相完全缺失Lo相（绿色区域）缺失Ld/Lo相分离（黄色区域），这是其最致命的缺陷 Martini 3的成功（K）：成功重现了Ld相区域（红色，左下角）成功重现了Lβ相区域（粉色/浅紫色，右下角）首次重现了Lo相区域（绿色，右上角高CHOL区域）成功捕捉到Ld/Lo共存带（黄色/橙色）成功捕捉到Lo/Lβ共存带（蓝色/青色）成功捕捉到Ld/Lβ共存带（深绿色）与实验相图的相似度达到定性一致，只在边界细节上有细微差异文章将模拟得到的相图与实验测定的三元相图逐点比较，验证了在不同DPPC/DOPC/CHOL组分比例下，Martini 3能够准确预测Ld、Lo、Lβ单相区以及它们的共存区，甚至捕捉到三相共存（Ld/Lo/Lβ）现象。这一成果证明了新的Martini 3脂质模型在捕捉复杂膜相行为方面的巨大进步，终于解决了粗粒化力场长达十余年无法准确描述脂质相分离的核心难题。 5. 模拟复杂生物膜与非层状结构真实细胞膜模型：作者使用新脂质组学构建了一个包含8种脂质、非对称分布的哺乳动物细胞质膜模型。该模型包含了胆固醇和鞘磷脂(SSM)等重要组分。模拟结果在膜厚、有序度、胆固醇翻转速率等方面都与Martini 2和全原子模拟的结果相符或更优。图3：Martini 3的复杂膜模拟。 (A,B) M3和C36模拟的系统快照。(C) M3、M2和C36模拟的各组分密度分布图。蛋白质-脂质相互作用：通过模拟钾离子通道Kir2.2和ADP/ATP载体等蛋白，证明了新模型能够准确识别蛋白质与特定脂质（如$PIP_2$和心磷脂）的结合位点。非层状相：新模型成功地模拟了DOPE脂质从层状到反相六方相 ($H_{II}$) 的转变，以及单油酸甘油酯 (MO) 自组装形成立方相 ($Q_{II}^D$) 的过程（图S6）。这些非层状结构在生物体内的膜融合过程以及作为药物递送载体（如脂质纳米粒, LNP）方面都至关重要。神经酰胺(Ceramide)和脂肪酸(Fatty Acid)的适用性：本文的框架为构建皮肤角质层脂质模型提供了坚实的基础。补充信息的全原子参考模拟中包含了神经酰胺(PCER, d18:1/16:0)的本体模拟，这为后续参数化提供了数据基础。同时，自动化脚本和灵活的映射方案使得构建不同链长的游离脂肪酸模型变得简单直接。更重要的是，通过精确重现胆固醇与磷脂的相分离行为，该工作验证了Martini 3中胆固醇模型的可靠性，这对于模拟由CER/CHOL/FFA组成的三元皮肤脂质体系至关重要。 Q&A Q1: 为什么新的映射方案能够区分仅相差2个碳原子的脂质链如此重要？ A1: 这个看似微小的改进是实现准确相行为预测的基石。原因如下：物理性质的差异：脂质尾链的长度直接决定了分子间的范德华相互作用强度和分子的几何形状。即使只相差两个碳原子（如DPPC的16碳链和DSPC的18碳链），也会导致它们的相变温度、膜厚度和堆积紧密程度产生显著差异。相分离的基础：在三元混合物中，胆固醇倾向于与更长、更直的饱和脂质链（如DPPC）紧密堆积形成有序的Lo相，而与带有扭结的不饱和脂质链（如DOPC）的相互作用较弱，后者形成无序的Ld相。如果模型无法从根本上区分不同长度的饱和链，就无法准确描述这种选择性的相互作用，也就无法重现相分离现象。化学特异性：能够区分细微的化学差异，是粗粒化模型从一个“通用”模型迈向“高保真”模型的关键一步，使其能够对更具体的生物化学问题做出可靠的预测。 Q2: 作者在参数化过程中提到了“人工在环优化 (human-in-the-loop)”，这具体是指什么？为什么不能完全自动化？ A2: “人工在环优化”是指在参数优化的过程中，研究人员需要根据多方面的、有时甚至是相互矛盾的验证结果，凭借专业知识和经验做出权衡与决策。在本文中，这意味着：多目标权衡：一个参数的改变可能会改善某个性质（如相变温度），但同时会恶化另一个性质（如单脂质面积APL）。例如，作者提到降低饱和尾链的角度力常数可以改善APL，但会导致$T_m$降低和相分离变差。自动化算法很难在这种多目标冲突中做出“科学上合理”的权衡。计算成本高昂：验证相分离或相变温度需要进行长时间的（数个微秒）模拟。将这样昂贵的计算嵌入一个全自动的优化循环（如贝叶斯优化）在计算上是不可行的。 “化学直觉”的引入：研究人员可以根据他们对物理化学原理的理解，有针对性地调整某些参数（如某个珠子的极性），而自动化算法通常是在整个参数空间中进行“黑箱”搜索，效率较低。 Q3: 新的Martini 3脂质组学如此成功，是否意味着全原子模拟不再重要了？ A3: 恰恰相反，这项工作更加凸显了全原子模拟的重要性。本文的成功是建立在一个多尺度的哲学之上的：全原子模拟是“老师”：Martini 3的参数化第一阶段，就是通过拟合高精度的CHARMM36全原子模拟数据来确定的。没有准确的全原子模拟作为“基准”，粗粒化模型的开发就成了无源之水。互补的角色：全原子模拟擅长提供精确的局部结构、相互作用能和短时动力学信息；而粗粒化模拟则擅长探索由这些局部相互作用涌现出的大尺度、长时间现象（如相分离）。两者是互补的，而非替代关系。未来的趋势是更多地将两者结合在多尺度工作流中。 Q4: 这项工作对于药物研发，特别是像mRNA疫苗这样的脂质纳米粒(LNP)递送系统，有什么意义？ A4: 意义非常重大。LNP的效率和稳定性与其内部的纳米结构密切相关，而这些结构往往是复杂的非层状相（如反相六方相或立方相）。本文展示了新的Martini 3模型能够准确模拟这些非层状相的形成。这意味着：配方筛选与优化：研究人员可以在计算机上高效地模拟由不同离子化脂质、辅助脂质和胆固醇组成的LNP配方，预测其内部结构，从而筛选出最有可能稳定包裹mRNA并高效递送的配方，大大缩短研发周期。机理研究：可以模拟LNP在不同生理环境（如内涵体的酸性环境）中的结构转变过程，从而在分子水平上理解其”内涵体逃逸”的关键机制。安全性评估：可以模拟LNP与细胞膜的相互作用，预测其潜在的细胞毒性或脱靶效应。关键结论与批判性总结潜在影响开启了大规模计算脂质组学：通过提供数千个经过验证的脂质模型和自动化工具，该工作将使广大研究人员能够以前所未有的规模和化学多样性来模拟复杂生物膜，从而推动“计算细胞生物学”的发展。解决了CG模拟的核心难题：成功地重现了脂质的相变和三元相分离，解决了长期困扰Martini力场的一个核心问题，极大地提升了其在研究膜微区、脂筏等生物学重要现象时的可靠性和预测能力。加速工业应用：通过提供能够模拟非层状相和复杂配方的工具，该工作将直接加速在药物递送（如LNP疫苗）、食品科学（如乳液稳定）和化妆品等领域的工业研发。研究局限性熵-焓补偿问题作为所有粗粒化模型的固有局限性，Martini 3仍然存在熵-焓补偿问题。这意味着其对温度的依赖性可能不完全准确，在远离参数化温度点（通常是310 K或323 K）时需谨慎使用。粗粒化过程中”自由度的减少”会导致焓和熵之间的平衡关系与全原子模拟不同，因此体系的热力学性质在较宽温度范围内的准确性有限。孔道形成能垒显著偏高这是Martini 3（以及所有当前Martini版本）面临的最严重的局限性之一。尽管在相行为描述上有显著改进，Martini 3模拟的膜上成孔的自由能垒仍然比全原子模拟高出数倍，这对研究涉及膜破坏的生物物理过程构成了重大障碍。定量证据：在补充信息图S18中，作者对比了Martini 3与全原子CHARMM36模拟DPPC双层膜成孔的自由能曲线： Martini 3计算的成孔自由能垒：约 170-180 kJ/mol CHARMM36全原子模拟的能垒：约 60-70 kJ/mol 差异：Martini 3的能垒几乎是全原子模拟的 2.5-3倍这一显著差异早在Bennett & Tieleman (2011) 的研究中就已被报道，当时对Martini 2和CHARMM36进行对比时发现了类似的问题。遗憾的是，即使经过Martini 3的全面改进，这一基本问题仍未得到解决。根本原因分析：这一问题的根源在于Martini力场对磷脂头基区域的简化表示： Q5珠子的化学非特异性：Martini使用单一的Q5型珠子来代表磷酸基团，这种高度简化的表示无法捕捉磷酸基团与水分子之间复杂的氢键网络和精细的静电相互作用缺失关键物理化学细节：成孔过程涉及磷脂头基的重新取向、水分子向疏水核心的渗透以及脂质尾链的复杂重排。这些过程对头基-水界面的精确描述极为敏感，而粗粒化模型在这方面天然存在局限熵效应的过度简化：成孔过程中的熵变（特别是水分子进入孔道时的构象熵和取向熵）在粗粒化模型中被显著低估对研究应用的影响：这一局限性使得Martini 3在以下研究场景中需要特别谨慎或不适用：电穿孔 (Electroporation)：在强电场下膜的击穿和孔道形成是该技术的核心，但能垒的严重高估会导致成孔时间尺度和阈值电场强度的预测完全偏离实际抗菌肽的膜破坏机制：许多抗菌肽通过形成跨膜孔道来杀死细菌，Martini 3可能无法正确捕捉这一过程的动力学和能量学膜融合的初期阶段：融合孔的形成和扩张是膜融合的关键步骤，能垒的偏差会影响对融合机制的理解去垢剂/表面活性剂的膜溶解：这类分子通过诱导膜缺陷和孔道来破坏脂质双层膜，Martini 3可能低估其效率未来改进方向：解决这一问题可能需要对磷酸基团及其周围水化层进行更精细的粗粒化处理，例如引入方向性相互作用或局部精细化策略。单脂质面积的系统性低估为了改善相变温度和相分离行为的预测，作者有意地将PC、PG和SM脂质的单脂质面积 (APL) 低估了约3 Å²。虽然这种”牺牲局部准确性以换取全局性质正确性”的策略在实践中是合理的，但它也意味着在研究对APL高度敏感的现象（如膜蛋白的镶嵌、膜张力的定量计算）时需要格外注意。蛋白质力场的兼容性虽然初步测试了与蛋白质的相互作用，但随着未来Martini 3蛋白质力场的进一步发展，脂质-蛋白质之间的相互作用参数可能需要重新评估和微调。目前的测试主要集中在已知的特异性结合（如$PIP_2$与离子通道），对于更复杂的蛋白质-膜相互作用（如膜曲率感应、蛋白质诱导的相分离）还需要更多验证。未来方向进一步扩大脂质库：将参数化范围扩展到更复杂的脂质，如糖脂、支链脂质和重要的信号脂质。自动化参数化：利用机器学习和自动化优化工具（如Swarm-CG）来进一步加速和完善新脂质的参数化流程，减少“人工在环”的需求。改进温度依赖性：探索开发具有温度依赖性势函数的新模型，以克服熵-焓补偿的限制，使其在更宽的温度范围内保持准确。小编笔记：对具体lipid类型，如ceramide，free fatty acid啥都没说学了几个新的指标，很好。有没有可能做一个Benchmark study，关于SC lipid的phase diagram以及和实验对？成孔自由能垒太高，这可咋办呀…做个新的工作来diss martini他们，甚至于调参来解决这个问题？

Molecular Dynamics · 2025-11-02

重塑细胞膜的关键角色：Martini 3粗粒化力场下的新一代胆固醇模型

重塑细胞膜的关键角色：Martini 3粗粒化力场下的新一代胆固醇模型本文信息标题: 用于胆固醇的Martini 3粗粒化力场作者: Luís Borges-Araújo, Ana C. Borges-Araújo, Tugba Nur Ozturk, Daniel P. Ramirez-Echemendia, Balázs Fábián, Timothy S. Carpenter, Sebastian Thallmair, Jonathan Barnoud, Helgi I. Ingólsson, Gerhard Hummer, D. Peter Tieleman, Siewert J. Marrink, Paulo C. T. Souza, and Manuel N. Melo 发表时间: 2023年10月5日单位: 里斯本新大学（葡萄牙），里昂大学（法国），劳伦斯利弗莫尔国家实验室（美国），卡尔加里大学（加拿大），马克斯·普朗克生物物理研究所（德国）等多个机构引用格式: Borges-Araújo, L., Borges-Araújo, A. C., Ozturk, T. N., Ramirez-Echemendia, D. P., Fábián, B., Carpenter, T. S., Thallmair, S., Barnoud, J., Ingólfsson, H. I., Hummer, G., Tieleman, D. P., Marrink, S. J., Souza, P. C. T., & Melo, M. N. (2023). Martini 3 Coarse-Grained Force Field for Cholesterol. Journal of Chemical Theory and Computation, 19(21), 7387–7404. https://doi.org/10.1021/acs.jctc.3c00547 摘要胆固醇通过调节脂质双层的流动性、刚性、通透性和组织结构，在生物膜中扮演着至关重要的角色。最新版本的Martini模型，即Martini 3，在相互作用平衡、分子堆积以及引入新型粒子类型和尺寸方面取得了显著改进。然而，新模型的发布也带来了对许多核心分子（包括胆固醇）进行重新参数化的需求。本文中，我们描述了一个Martini 3胆固醇模型的开发和验证过程，解决了与其键合设置、形状、体积和疏水性相关的问题。我们提出的新模型缓解了其Martini 2前身的一些局限性，同时保持或改善了其整体行为。核心结论成功开发并验证了一款新的Martini 3胆固醇粗粒化模型。该模型在形状、疏水性和动力学稳定性方面均有显著提升。通过创新的“单框架虚拟位点”拓扑结构，彻底解决了Martini 2模型中存在的、由LINCS约束算法导致的“人工温度梯度”artifact。新模型更准确地再现了胆固醇的物理化学性质。它修正了旧模型过于“粘稠”（过度亲脂）的问题，其形状也更逼真，从而在模拟中实现了更准确的分子堆积。新模型在多种复杂生物场景中表现优异。它能正确再现胆固醇对膜的“增稠”和“致密”效应、在三元脂质体系中诱导相分离，并能准确识别其在多个重要膜蛋白上的结合位点。背景胆固醇是动物细胞膜中不可或缺的“万能调解员”。它像楔子一样插入磷脂分子之间，灵巧地调节着细胞膜的流动性、刚性和通透性。同时，它还是形成“脂筏”——一种富含特定脂质和蛋白质的微观区域——的关键驱动力，深刻影响着细胞信号转导等多种生命过程。此外，胆固醇还能直接与膜蛋白相互作用，调控其功能，并且是合成类固醇激素的前体。近年来，随着mRNA疫苗等基因疗法的发展，胆固醇作为脂质纳米颗粒递送系统的关键组分，其重要性愈发凸显。为了在原子尺度下理解这些复杂的生物物理过程，分子动力学 (MD) 模拟已成为不可或缺的研究工具。然而，全原子模拟的计算成本极高，难以企及细胞膜重塑、相分离等发生在大尺度（微米级）和长时程（毫秒级）上的现象。为此，粗粒化 (Coarse-Grained, CG) 模型应运而生。其中，Martini力场将约4个重原子简化为一个CG粒子（bead），极大地提升了模拟效率，已成为生物膜模拟领域最流行的CG力场之一。然而，广泛使用的Martini 2版本存在一些固有缺陷。特别是对于蛋白质和胆固醇这类环状刚性分子，模型会表现出过度的疏水性和自相互作用，即过于“粘稠”。此外，Martini 2的胆固醇模型在使用GROMACS中的LINCS约束算法时，会产生人工的温度梯度，即不同分子（如胆固醇和磷脂）在模拟中会表现出不同的温度，这是一个严重的物理artifact。2021年发布的全新Martini 3框架通过引入更多样的粒子类型和尺寸，并优化相互作用平衡，系统性地解决了这些问题。但这也意味着，包括胆固醇在内的几乎所有分子都需要重新进行参数化。关键科学问题本研究的核心科学问题是：如何构建一个全新的、与Martini 3框架兼容的胆固醇粗粒化模型，该模型不仅能解决Martini 2版本中存在的数值不稳定（温度artifact）和物理不准确（过度疏水）等关键问题，还能在更广泛的生物物理场景中准确地再现胆固醇的结构、热力学和动力学行为？具体来说，研究团队需要攻克以下几个技术难点：拓扑结构设计：如何设计一个既能精确描述胆固醇刚性环状结构，又能在数值上保持稳定、与常用约束算法（如LINCS和CCMA）良好兼容的键合网络？化学性质校准：如何通过精细地选择CG粒子类型，来修正胆固醇的疏水性，使其在水/油两相中的分配行为与实验值相符？形状与堆积：如何让简化的CG模型能够再现胆固醇独特的、带有“粗糙”面（有甲基伸出）和“光滑”面的三维形状，从而实现其在脂质膜中正确的堆积和组织行为？综合性能验证：新模型能否在多种复杂的膜环境中（不同饱和度的脂质、三元混合物相分离、与蛋白质相互作用等）都表现出优于或至少不逊于旧模型的性能？创新点创新的单框架虚拟位点拓扑：设计了一种新颖的“单框架虚拟位点 (single-frame virtual site)”拓扑结构。这一设计巧妙地解决了Martini 2模型中因“双框架”结构与LINCS约束算法不兼容而产生的人工温度梯度artifact，同时保证了模型在长时程模拟中的稳定性。更逼真的分子形状与化学性质：通过引入新的“微小 (tiny)”尺寸粒子来显式地表示胆固醇环上的两个轴向甲基，并精心组合不同类型的CG粒子，新模型在三维形状（如溶剂可及表面积）和疏水性（如油水分配自由能）上都更接近全原子参考和实验值。跨平台兼容性：新的拓扑结构不仅解决了GROMACS中的LINCS问题，还天然兼容OpenMM模拟引擎中的CCMA约束算法，而后者无法稳定模拟Martini 2的胆固醇模型。这极大地增强了新模型在不同计算化学社区中的通用性。全面而严苛的验证：新模型经历了一场“全能大考”，系统性地验证了其在再现胆固醇诱导的膜增厚、面积压缩、脂质排序、在复杂三元体系中的相分离行为，以及与三种不同类型膜蛋白（GPCRs和离子通道）的相互作用等多种关键生物物理现象中的表现，证明了其广泛的适用性和可靠性。研究内容核心方法论：构建新一代Martini 3胆固醇模型构建一个优秀的粗粒化模型，如同创作一幅神似的写意画，既要抓住精髓，又要舍弃繁琐。作者采用了一套自下而上、反复迭代的参数化流程，每一步都以高精度的全原子模拟数据为“金标准”。 graph TD subgraph "Martini 3 胆固醇模型参数化流程" direction LR A("1.建立参考体系 长时间全原子模拟 (CHARMM36力场, >1µs)"); A --> B["2.CG映射与拓扑设计 确定粒子数量、位置和连接方式 (创新的'单框架虚拟位点')"]; B --> C["3.优化键合参数 匹配键长、键角、二面角分布 (对比CG与AA映射后的分布)"]; C --> D["4.优化非键参数 (粒子类型选择) 匹配热力学性质 (如油水分配自由能)"]; D --> E{"5.综合性能验证 (膜性质、相分离、蛋白相互作用等)"}; E -- "不满足要求" --> B; E -- "满足要求" --> F("最终模型"); end 1. 模拟设置与分析工具参考标准：所有粗粒化模型的开发都以CHARMM36全原子 (AA) 力场的模拟结果为基准。AA模拟的时长至少为1微秒，以确保充分的采样。粗粒化模拟：CG模拟使用GROMACS或OpenMM进行，时长通常在10微秒以上，以检验模型的长期稳定性和物理行为。分析软件：整个流程广泛使用了多种Python科学计算库，如MDAnalysis用于轨迹分析，Voro++用于计算单位脂质面积，pymbar用于自由能计算，LiPyphilic和PyLipID分别用于分析胆固醇翻转和停留时间。 2. 更逼真的映射与形状图1：Martini 3胆固醇模型的参数化。(a) 化学结构与粗粒化映射。(b) 新模型的Connolly表面与全原子参考对比。(c) Martini 2（双框架）与Martini 3（单框架）虚拟位点拓扑示意图。(d) 溶剂可及表面积(SASA)对比。(e) 均方根偏差(RMSD)对比。(f, g) 辛醇/水和十六烷/水分配自由能对比。映射方案：新模型将胆固醇分子简化为9个CG粒子。例如，根据附录中的Table S4，代表柔性尾链末端的C2粒子，实际上对应着全原子模型中的C23, C24, C25, C26, C27共5个碳原子。形状优化：为了更准确地描述胆固醇独特的、带有“粗糙”面（有甲基伸出）和“光滑”面的三维形状，作者创新地使用了两个**“微小 (tiny)”**尺寸的粒子 (R5, R6) 来显式地表示从甾环平面伸出的两个轴向甲基。这使得新模型的溶剂可及表面积 (SASA) 和整体形状都与全原子参考更为接近。 3. 解决数值稳定性的“单框架”拓扑 Martini 2的问题：旧模型使用“双框架虚拟位点”拓扑来维持刚性。它由两个共享一条边的三角形约束框架构成，像一个可以折叠的铰链。这种设计在GROMACS的LINCS约束算法下容易出现收敛问题，导致能量无法在分子内正确传递，从而产生胆固醇分子“过冷”的人工温度梯度artifact。 Martini 3的解决方案：新模型采用更简洁的“单框架虚拟位点”拓扑。它仅使用R1, R2, C1三个粒子构成一个单一的刚性三角形约束框架，其余的甾环粒子（R3, R4, R5, R6）则作为无质量的虚拟位点，其位置由这个框架的三个顶点唯一几何确定。为了保持质心不变，这些虚拟位点的质量被重新分配到了三个框架粒子上。图S2：Martini 3胆固醇模型的温度差异。附录中的这张图是关键证据，它显示了在一个包含DLIPC、DPPC和胆固醇的混合体系中，使用新模型模拟时，三种分子的平均温度（柱状图a）和瞬时温度（曲线图b）都稳定在设定的300K附近，完全消除了Martini 2模型中存在的温度梯度artifact。 4. 更平衡的化学性质修正过度疏水性：Martini 2胆固醇模型过于“粘稠”，其油水分配自由能远高于实验值。Martini 3模型通过精心组合不同化学性质的粒子类型来解决此问题：甾环核心 (R2, R3, R4) 使用SC3类型粒子，这类粒子被设计用于环烷烃，疏水性适中。伸出的甲基和烷基尾链 (R5, R6, C1, C2) 使用TC2和C2类型粒子，它们被设计用于支链烷烃，与饱和脂质尾链（C1类型）有良好的相互作用。验证结果：通过自由能微扰方法计算，新模型的辛醇/水和十六烷/水分配自由能与实验或理论参考值的吻合度都得到了显著提升。结果与分析：新模型的全面性能验证新模型在一系列严苛的测试中展现了其优越的性能，证明了其在多种生物物理场景下的可靠性。 1. 在脂质膜中的基本行为图2：胆固醇在不同脂质双层中的插入行为。(a) 胆固醇羟基(ROH)的密度分布图。(b) 胆固醇在不同饱和度脂质膜中的跨膜翻转（flip-flop）速率。正确的膜内定位与翻转：在饱和脂质膜（如DPPC）中，新模型能像真实胆固醇一样，以经典的“直立”姿态插入膜中，羟基锚定在磷酸头基区域。随着膜不饱和度的增加，模型也开始出现少量平行于膜中心排列的非标准构象，并表现出翻转速率随不饱和度增加而加快的趋势，这与实验观察和物理预期一致。 2. 对膜物理性质的调控作用图4：胆固醇对DPPC和POPC双层膜的影响。(a, d) 膜厚度变化。(b, e) 单位脂质面积(APL)变化。(c, f) 脂质尾链平均有序度(S-order)变化。经典的“增稠”与“致密”效应：与实验和全原子模拟一致，随着胆固醇浓度的增加，新模型能够正确地使DPPC（饱和）和POPC（不饱和）膜增厚，同时压缩脂质分子，使其平均占据的面积(APL)减小。强大的“排序”能力：胆固醇的加入显著增加了脂质尾链的有序度（S-order），即让原本杂乱的尾链变得更加挺直有序。S-order的计算公式为： $S = \frac{1}{2}(3\langle(\cos\theta)^2\rangle - 1)$ 其中θ是CG粒子对之间的连线与膜法线（z轴）的夹角。新模型能很好地再现这一排序效应。跨平台一致性：附录中的图S8显示，使用GROMACS和OpenMM两种不同的模拟软件，新模型在预测膜厚度、APL和有序度等性质时，给出了几乎完全一致的结果，这强有力地证明了新模型的跨平台兼容性和稳健性。 3. 诱导相分离的能力图5：胆固醇对三元脂质体系相分离的影响。比较了Martini 2 (a-d) 和Martini 3 (e-h) 模型在不同温度下诱导相分离的能力。(i, j) 定量分析了DPPC-DPPC和CHOL-DPPC的接触分数，分数越高表示分离越明显。再现液有序相：在由饱和脂质(DPPC)、不饱和脂质(DLIPC)和胆固醇构成的三元体系中，新模型成功地再现了相分离现象：胆固醇倾向于与DPPC聚集，形成致密的液有序(Lo)相，而DLIPC则形成液无序(Ld)相。优于旧模型：定量分析显示，在不依赖温度artifact的情况下，新模型诱导相分离的能力与Martini 2相当甚至略有改善。虽然对于更难分离的DPPC/DOPC/CHOL体系，新旧模型都表现不佳（这被归因于脂质模型本身的问题），但新模型至少为研究复杂的细胞膜组织行为提供了一个更可靠的出发点。 4. 与膜蛋白的相互作用研究者进一步测试了新模型与三种重要的膜蛋白（β2肾上腺素受体、SMO受体和VDAC1离子通道）的相互作用。图6：胆固醇与β2AR的识别和结合。 (a) 胆固醇的占据密度图。(b) 蛋白表面按胆固醇停留时间着色。(c) 实验晶体结构中发现的胆固醇。(d) 模拟快照显示胆固醇结合在已知位点。图7：胆固醇与SMO的识别和结合。精准识别结合位点：在长时间的模拟中，新模型能够准确地识别并稳定结合到这些蛋白上已知的、由实验（如X射线晶体学）或全原子模拟确定的胆固醇结合位点上。更真实的动力学：相比Martini 2模型由于过度粘稠而导致的微秒级停留时间，新模型的胆固醇与蛋白的相互作用更加动态，停留时间在纳秒级，虽然可能略有低估，但通过快速的交换，依然能维持在高占据率的结合位点上。这为研究胆固醇对膜蛋白功能的动态调控提供了更真实的视角。 Q\&A Q1: Martini 2的胆固醇模型有什么具体问题？Martini 3是如何从根本上解决的？ A1: Martini 2模型主要有两个核心问题： 1. 数值不稳定性（温度artifact）：它采用的“双框架虚拟位点”拓扑结构，在GROMACS的LINCS约束算法下容易出现收敛失败。这导致动能无法在分子内部正确分配，使得胆固醇分子自身的温度会显著低于体系的设定温度，这是一个严重的物理artifact，会影响相分离等性质。Martini 3通过设计更简洁、更稳固的**“单框架虚拟位点”拓扑**，从根本上解决了这个问题。 2. 物理不准确性（过度疏水）：Martini 2的粒子类型和相互作用定义使得胆固醇分子过于“粘稠”，即它与疏水环境（如脂质尾链）的相互作用过强，而与水相的排斥也过强。这导致其油水分配自由能与实验值偏差很大。Martini 3通过引入更多样化的粒子类型（如SC3, TC2, C2）并重新优化它们的组合，更精细地刻画了胆固醇不同部分的化学性质，使其整体疏水性回归到更合理的水平。 Q2: 什么是“虚拟位点 (Virtual Site)”，为什么在胆固醇这类刚性分子的粗粒化模型中要使用它？ A2: “虚拟位点”是一个在模拟中没有质量的粒子，它的坐标不是通过积分运动方程得到的，而是根据体系中其他“真实”粒子的位置实时计算出来的。在粗粒化胆固醇模型中使用虚拟位点主要有两个目的： 1. 维持刚性结构：胆固醇的甾环是一个非常刚性的结构。如果用大量的键和角来维持这个形状，会导致模型中出现高频振动，迫使模拟使用很小的时间步长，从而失去粗粒化的速度优势。通过定义一个由少数真实粒子构成的刚性框架（如“单框架”中的三角形），然后将其他粒子定义为基于这个框架计算出的虚拟位点，就可以在保持整体刚性的同时，避免引入过多的键合相互作用。 2. 提高数值稳定性：如前所述，一个设计良好的虚拟位点拓扑结构可以避免与约束算法的冲突，提高模拟的稳定性和准确性。 Q3: 新模型在膜相分离的模拟中表现如何？有什么改进和仍然存在的挑战？ A3: 新模型在相分离方面的表现可以说是有显著进步，但仍有提升空间。进步之处：它成功地再现了DPPC/DLIPC/CHOL三元体系的相分离。更重要的是，它是在没有温度artifact的情况下实现这一点的。而Martini 2模型之所以能看到相分离，部分原因是由于胆固醇“过冷”这一artifact增强了其与DPPC的聚集。因此，Martini 3的成功是基于更正确的物理基础。挑战之处：对于更难分离的DPPC/DOPC/CHOL体系，新模型未能观察到预期的相分离。但作者指出，这很可能不是胆固醇模型本身的问题，而是因为当前Martini 3的DPPC/DOPC脂质对模型本身就难以相分离。这说明，一个体系的准确模拟依赖于力场中所有组分的共同努力，对胆固醇的改进还需要未来对脂质模型的进一步优化来相辅相成。 Q4: 论文提到新模型在OpenMM中也能稳定运行，这有什么重要意义？ A4: 这一点具有非常重要的实践意义。不同的MD模拟引擎使用不同的算法来处理键合约束。例如，GROMACS主要使用LINCS算法，而OpenMM则常用CCMA算法。Martini 2胆固醇模型的“双框架”拓扑与CCMA算法不兼容，导致其在OpenMM中无法稳定运行。而Martini 3胆固醇模型采用的“单框架”设计，既解决了GROMACS中的LINCS问题，又天然兼容OpenMM的CCMA算法，如附录图S8所示，两种软件给出的结果几乎完全一致。这极大地增强了模型的可用性和在不同科研社区间的通用性。关键结论与批判性总结潜在影响提升模拟可靠性：通过解决关键的技术artifact并提高物理准确性，这款新的Martini 3胆固醇模型为整个生物膜模拟领域提供了一个更可靠、更稳健的基础工具，将提升大量依赖于该模型的下游研究（如脂筏、病毒包膜、脂质纳米颗粒等）的质量。促进多平台协作：解决了跨主流MD引擎的兼容性问题，有助于统一不同实验室和研究社区的模拟标准，促进结果的可重复性和比较。加速药物研发：一个更准确的胆固醇模型对于模拟其与GPCRs等药物靶点的相互作用至关重要，有助于更精确地理解药物的变构调节机制和设计靶向特定脂质环境的药物。研究局限性部分性质仍有偏差：尽管取得了巨大进步，新模型在某些定量性质上仍非完美。例如，它仍然略微低估了胆固醇对膜的增厚效应，并且在高度不饱和的膜中，其跨膜翻转速率可能被高估。依赖于其他模型：胆固醇在膜中的行为（如相分离）强烈依赖于与之相互作用的脂质模型。当前模型在某些三元体系中的表现不佳，凸显了其性能受限于整个Martini 3脂质力场的整体发展水平。动力学校准的挑战：粗粒化模型由于表面光滑，动力学过程通常会比全原子模拟快4倍左右。虽然这是一个已知的特征，但对于需要精确动力学信息的场景（如计算解离速率），仍需谨慎使用或进行额外校准。未来方向力场的协同进化：未来的工作将集中于对Martini 3的磷脂模型进行迭代改进，以解决与胆固醇相互作用时表现出的剩余偏差（如相分离问题），实现整个脂质力场的协同优化。拓展到其他甾醇：利用本次参数化建立的成功经验和拓扑设计，可以将其推广到其他重要的甾醇分子，如植物甾醇、麦角固醇（真菌）和hopanoids（细菌），构建一个完整的Martini 3甾醇家族。更复杂的应用验证：将新模型应用于更具挑战性的生物系统中，例如模拟真实细胞器（如内质网）膜的复杂脂质组成、病毒与宿主细胞膜的融合过程，或包含多种膜蛋白和脂筏的拥挤细胞膜环境。

Molecular Dynamics · 2025-11-02

Martini 3 脂质组学补充材料概览：方法、验证与应用

Martini 3 脂质组学补充材料概览：方法、验证与应用本文档是对Martini 3脂质组学论文（Souza et al., 2021, JACS Au）补充材料的系统性总结。补充材料共61页，包含详细的验证实验、方法学说明及模型局限性讨论。补充结果概述 A. 双层膜弯曲模量的改进研究问题：Martini 2系列模型系统性地高估了脂质双层膜的弯曲模量（$k_c$），这影响了膜变形和膜重塑过程的模拟准确性。方法：实空间起伏法（RSF）：通过分析膜表面高度起伏的功率谱计算$k_c$ 屈曲法（Buckling）：对小尺寸膜片施加表面张力，通过屈曲转变计算$k_c$ 关键发现： Martini 3在弯曲模量精度上显著优于Martini 2，多数脂质的$k_c$值更接近实验数据 POPC的$k_c$从Martini 2的约40-50 $k_BT$降低至Martini 3的约20-30 $k_BT$（实验值约18-25 $k_BT$）不同计算方法（RSF vs Buckling）给出的结果基本一致，验证了参数化的稳健性物理意义：更准确的弯曲模量使得Martini 3能够更好地模拟膜融合、内吞、出胞等生物学过程。 B. 自动生成脂质拓扑及双层膜性质探索研究目的：展示Martini 3的自动化工作流程，系统性地生成并验证大量脂质的拓扑参数。方法：使用自动化脚本从化学结构生成Martini 3脂质拓扑对每种脂质进行标准双层膜模拟（NPT系综，323 K）计算关键物理量：面积密度（APL）、双层厚度（$d_{HH}$）、序参数（$S_{CD}$）、相变温度（$T_m$）关键发现：成功生成并验证了数百种脂质分子的拓扑多数脂质的APL、厚度等性质与实验数据吻合良好发现了一些系统性偏差：某些长链饱和脂质的$T_m$略高于实验值工具化成果：这一自动化流程已集成到insane.py工具和Martini Lipidome Database中，用户可以快速构建含有任意脂质组成的膜体系。 C. 中性脂质的密度和界面张力研究对象：中性脂质（如二酰基甘油DAG、三酰基甘油TAG、胆固醇酯CE等）在膜结构和脂滴形成中起重要作用。验证指标：体密度：纯相中性脂质的密度界面张力：中性脂质与水的界面张力关键发现： Martini 3对中性脂质的密度再现良好，与实验值的偏差在5%以内界面张力的预测也较为准确，特别是TAG和CE的水-脂界面性质这些参数对于模拟脂滴形成、脂筏结构等现象至关重要应用前景：为研究脂质代谢、脂滴动力学提供了可靠的力场基础。 D. 离子调控的磷脂酰丝氨酸相分离生物学背景：磷脂酰丝氨酸（PS）是重要的阴离子脂质，其在细胞膜中的分布受到离子（特别是Ca²⁺）的调控。模拟设计：构建POPC/POPS混合膜体系改变溶液中Ca²⁺浓度观察PS的相分离行为关键发现：高浓度Ca²⁺能够诱导PS富集区域的形成（相分离） Martini 3能够再现PS-Ca²⁺的特异性相互作用相分离的程度与Ca²⁺浓度呈正相关生物学意义：PS的相分离与细胞信号转导、膜融合等过程密切相关，Martini 3为研究这些现象提供了工具。 E. 非层状脂质相模拟研究背景：某些脂质在特定条件下会形成非层状相，如反向六方相（HII）、立方相（QIID）等，这些相在膜融合和膜蛋白功能中有重要作用。验证体系：反向六方相（HII）：DOPE（二油酰基磷脂酰乙醇胺）立方相（QIID）：单油酸甘油酯（MOG）关键发现： Martini 3能够自发形成并稳定HII相，与实验观察一致立方相的形成也得到了初步验证非层状相的形成温度和相转变温度与实验数据基本吻合技术挑战：非层状相的模拟对体系尺寸和平衡时间要求较高，需要数微秒级别的模拟才能充分平衡。 F. 真实脂质组成的复杂膜模拟研究目的：验证Martini 3在生理相关的复杂膜体系中的表现。模拟体系：类质膜（plasma membrane-like）：包含PC、PE、PS、胆固醇等多种组分线粒体膜：富含心磷脂（cardiolipin）细菌膜：包含特殊脂质如脂多糖（LPS）关键发现： Martini 3能够稳定模拟包含10种以上不同脂质的复杂膜膜的整体厚度、流动性等性质与实验数据一致观察到了脂筏样结构（胆固醇富集区）的自发形成应用价值：为研究膜的横向组织、蛋白质的膜定位提供了更真实的环境。 G. 蛋白质-脂质相互作用研究研究问题：蛋白质如何影响膜的局部结构？Martini 3能否准确捕捉蛋白质-脂质的特异性相互作用？验证体系：跨膜蛋白：如GPCR、离子通道外周膜蛋白：如annexin、PH结构域关键发现： Martini 3能够再现蛋白质对膜厚度的扰动（hydrophobic mismatch效应）特定脂质（如PIP2）在蛋白质周围的富集现象得到了正确描述外周膜蛋白的膜结合取向与实验/全原子模拟一致技术要点：蛋白质使用Martinize2工具转换为粗粒化模型，保持与脂质力场的兼容性。模型局限性与未来方向（Supplementary Discussion H）当前局限性熵-焓补偿问题： Martini 3通过调整LJ势能参数来匹配实验观测，但这种做法可能导致熵和焓的贡献不完全正确例如，某些相变温度是通过调整相互作用强度得到的，而非通过正确的微观机制孔道形成能垒： Martini模型中膜的孔道形成自由能垒偏低，导致大分子（如DNA、蛋白质）更容易穿膜这可能影响膜通透性和跨膜传输过程的模拟电荷相互作用的处理：粗粒化模型中电荷的有效性需要进一步优化特别是在多价离子（如Ca²⁺、Mg²⁺）存在时，相互作用的精度有待提高特定脂质的参数化：一些特殊脂质（如含有不饱和键的脂质、含有糖基的糖脂等）的参数仍需进一步优化长链饱和脂质的相变温度系统性偏高未来改进方向开发更精细的粗粒化策略（如超粗粒化、多尺度耦合）引入极化效应以更准确描述电荷相互作用针对特定生物学问题（如膜融合、内吞）进行专门的参数优化与实验（特别是中子散射、冷冻电镜）更紧密结合，提供更多验证数据方法学要点（Supplementary Methods I-M） I. 参考模拟（Reference Simulations）目的：建立标准化的模拟协议，确保不同研究者能够复现结果。标准流程：体系构建：使用insane.py生成初始结构能量最小化：最速下降法，$F_{max} < 10$ kJ·mol⁻¹·nm⁻¹ 平衡模拟：NVT（100 ps）→ NPT（1 ns），逐步释放位置约束生产模拟：NPT系综，半各向同性压力耦合，时间步长20 fs 关键参数：温度：323 K（v-rescale恒温器，τ_T = 1.0 ps）压力：1 bar（Parrinello-Rahman压力耦合，τ_P = 12.0 ps）静电：反应场（RF），截断1.1 nm 范德华：势能平移（potential-shift），截断1.1 nm J. 实验基准验证（MIB - Martini lipid Benchmark） MIB数据库：系统性收集了文献中报道的脂质双层膜实验数据，包括：面积密度（APL）双层厚度（$d_{HH}$）序参数（$S_{CD}$）相变温度（$T_m$）验证流程：对每种脂质进行标准模拟计算上述物理量与MIB数据库中的实验值对比量化模型的系统性偏差统计指标：平均绝对误差（MAE）均方根误差（RMSE） Pearson相关系数 K. 复杂双层膜的构建工具：insane.py脚本支持的功能：任意脂质组成：可指定每种脂质的比例不对称膜：上下叶片可以有不同的脂质组成嵌入蛋白质：自动在膜中插入粗粒化蛋白质溶剂离子：自动添加水和盐使用示例： insane.py -l POPC:70 -l CHOL:30 -prot protein.pdb -sol W -salt 0.15 -o system.gro L. 蛋白质-脂质相互作用的建模蛋白质粗粒化：使用Martinize2工具将全原子蛋白质结构转换为Martini模型保持二级结构稳定（通过弹性网络或Go模型）膜嵌入：根据蛋白质的疏水性残基分布确定跨膜区域使用insane.py自动将蛋白质嵌入膜中模拟策略：初始阶段对蛋白质施加位置约束，让脂质充分弛豫逐步释放约束，观察蛋白质-脂质的动态相互作用 M. 相行为的定量分析 Lindemann指数：用于判断脂质尾链的有序-无序转变（凝胶相-流体相） \[\delta_i = \frac{1}{N-1} \sum_{j \neq i} \frac{\sqrt{\langle r_{ij}^2 \rangle - \langle r_{ij} \rangle^2}}{\langle r_{ij} \rangle}\] $\delta_i < 0.1$：有序相（凝胶相） $\delta_i > 0.1$：无序相（流体相）富集指数：用于定量描述脂质相分离程度 \[E_A = \frac{N_A^{local} / N_{total}^{local}}{N_A^{global} / N_{total}^{global}}\] $E_A > 1$：脂质A在局部富集 $E_A < 1$：脂质A在局部贫化数据资源（Supplementary Data N） Martini Lipidome Database 内容： 500+ 脂质分子的拓扑文件（.itp格式）每种脂质的验证数据（APL、厚度、相变温度等）标准化的命名规则和分类系统访问方式：在线数据库：cgmartini.nl/lipidome GitHub仓库：包含所有拓扑文件和示例脚本 API接口：提供Python API，方便自动化工作流程支持批量下载和参数查询应用示例： from martini_lipidome import Lipid # 获取POPC的拓扑信息 popc = Lipid('POPC') print(popc.area_per_lipid) # 输出：0.61 nm² print(popc.phase_transition_temp) # 输出：271 K 总结本补充材料为Martini 3脂质组学的开发和验证提供了全面、系统的技术文档。关键要点包括：方法学创新：两阶段参数化策略（阶段I：单体性质，阶段II：凝聚相性质）确保了模型的物理合理性大规模验证：通过MIB基准数据库对数百种脂质进行了系统性验证，量化了模型的精度和局限性工具化与开放：提供了完整的工具链（insane.py、Martinize2、Lipidome Database）和API，降低了使用门槛应用导向：针对复杂膜体系、蛋白质-脂质相互作用等实际应用场景进行了专门优化透明的局限性讨论：明确指出了模型当前的不足（如熵-焓补偿、孔道形成能垒等），为未来改进指明了方向展望：Martini 3为膜生物学、药物递送、膜蛋白功能等研究提供了强大的模拟工具。随着参数的持续优化和新功能的开发（如极化模型、多尺度耦合），其应用范围将进一步扩大。参考文献 Souza, P. C. T.; Alessandri, R.; Barnoud, J.; Thallmair, S.; Faustino, I.; Grünewald, F.; Patmanidis, I.; Abdizadeh, H.; Bruininks, B. M. H.; Wassenaar, T. A.; Kroon, P. C.; Melcr, J.; Nieto, V.; Corradi, V.; Khan, H. M.; Domański, J.; Javanainen, M.; Martinez-Seara, H.; Reuter, N.; Best, R. B.; Vattulainen, I.; Monticelli, L.; Periole, X.; Tieleman, D. P.; de Vries, A. H.; Marrink, S. J. Martini 3: A General Purpose Force Field for Coarse-Grained Molecular Dynamics. JACS Au 2021, 1 (6), 587–608. https://doi.org/10.1021/jacsau.1c00203 文档说明：本文档基于Martini 3脂质组学论文的补充材料（oc5c00755_si_001.pdf，共61页）整理而成，旨在为读者提供快速、系统的技术概览。详细数据和图表请参考原始补充材料。

Molecular Dynamics · 2025-11-02

Martini 3珠子类型与命名规则：粗粒化分子动力学力场的完整指南

title: “Martini 3 Bead Types and Naming Conventions: A Comprehensive Guide” date: “2025-05-27” description: “Martini 3 珠子类型与命名规则的完整指南。详细介绍粗粒化分子动力学力场的珠子类型系统，包括命名规范、参数设置和应用建议。” tags: [martini3, coarse-grained, molecular-dynamics, force-field, bead-types, parametrization, cg-modeling] thumbnail: “/assets/img/thumbnail_mine/wh-m992d8.jpg” image: “/assets/img/thumbnail_mine/wh-m992d8.jpg” — 主要参考资料： https://doi.org/10.1038/s41592-021-01098-3 Supporting information for: Martini 3: A General Purpose Force Field for Coarse-Grained Molecular Dynamics https://github.com/ricalessandri/Martini3-small-molecules/blob/main/tutorials/building_block_table.pdf https://advanced.onlinelibrary.wiley.com/doi/full/10.1002/adts.202100391 https://cgmartini.nl/docs/tutorials/Martini3/Small_Molecule_Parametrization/ 1. 引言 (Introduction) Martini 力场是一种广泛应用于生物分子模拟的粗粒化 (Coarse-Grained, CG) 模型 (1)。近年来，经过大幅改进和重新参数化的 Martini 3 版本正式发布 (1)。Martini 3 旨在提供一个通用性更强的 CG 力场，不仅适用于脂质、蛋白质、核酸和糖类等生物大分子体系 (4)，也拓展到了对多种小分子、碳纳米材料以及聚合物的研究 (7)。相较于早期版本，Martini 3 的核心改进包括更优化的非键相互作用平衡、引入了新的珠子 (bead) 类型（包括不同尺寸和化学特性的珠子）、并增强了对特定相互作用（如氢键和电子极化效应）的描述能力 (1)。这些改进使得 Martini 3 能够更准确地预测分子的堆积模式和相互作用，从而在更广泛的应用领域中提供可靠的模拟结果 (1)。Martini 模型通常采用“四对一”的映射方案，即平均四个重原子及其相连的氢原子被粗粒化为一个相互作用中心（珠子），但对于环状结构等特殊化学基团，也支持更高分辨率的映射 (2)。本报告旨在详细阐述 Martini 3 力场中珠子的类型、命名方式的传统和原则，并深入探讨其参数化策略和分子映射方法。最后，将通过一个具体的聚合物——聚[2-(N-氧化-N,N-二乙基氨基)甲基丙烯酸乙酯] (poly[2-(N-oxide-N,N-diethylamino)ethyl methacrylate])——的映射实例，展示如何将这些理论知识应用于实践。 2. Martini 3 核心珠子 (Bead) 类型与命名传统 (Martini 3 Core Bead Types and Naming Conventions) Martini 3 模型的基石在于其多样化的珠子类型，这些珠子代表了不同化学性质的分子片段。理解这些珠子的分类、尺寸和命名规则对于正确构建和解读 CG 模型至关重要。 2.1 主要珠子类型 (Main Bead Types) 与早期版本类似，Martini 3 保留了基于极性的四种主要珠子类型 (8)： P (Polar): 极性珠子，代表强极性基团。 N (Non-polar/Intermediate polarity): 中等极性或非极性珠子，代表具有一定极性或非极性的基团。 C (Apolar/Carbon-like): 非极性珠子，通常代表疏水性的烷烃链等。 Q (Charged): 带电荷珠子，代表离子化的基团。在 Martini 3 中，这些主要类型得到了扩展和细化，引入了新的专用珠子类型 (8)： W (Water): 特定的水珠子，与 Martini 2 中的极性 P4 珠子不同，W 珠子经过独立优化，避免了旧模型中水在室温下结冰等问题。 D (Divalent ions): 二价离子珠子。 X (Halo-compounds): 用于描述含卤素原子的基团。这些主要类型（P, N, C, Q, X）进一步划分为多个亚型，通过数字后缀（通常为1-6）表示其相对极性程度或相互作用强度，数字越大通常表示极性越强或相互作用越强 (10)。Martini 3 将可能的相互作用能级从 Martini 2 的10个扩展到了22个，从而可以更精细地描述不同化学基团间的相互作用差异 (8)。此外，还引入了字母后缀来表征特定的化学特性，如氢键给体/受体能力和电子极化效应 (8)。 2.2 珠子尺寸 (Bead Sizes) Martini 3 引入了三种不同尺寸的珠子，以适应不同分辨率的粗粒化需求，这对于精确描述分子形状和堆积至关重要 (8)： Regular (R): 常规尺寸珠子，其 Lennard-Jones (LJ) 参数中的$\sigma$值约为 0.47 nm。设计用于标准的“4对1”映射，即代表4个重原子及其相连的氢原子。 Small (S): 小尺寸珠子，$\sigma$值约为 0.41 nm。设计用于“3对1”的映射，即代表3个重原子。 Tiny (T): 微小尺寸珠子，$\sigma$值约为 0.34 nm。设计用于“2对1”的映射，即代表2个重原子。这三种尺寸的珠子之间的交叉相互作用 (R-S, R-T, S-T) 都经过了专门的参数化，以确保整个力场的平衡性 (8)。小尺寸 (S) 和微小尺寸 (T) 珠子特别适用于描述环状结构（如芳香环和脂肪环）以及其他需要更高分辨率的线性或支链化学基团 (4)。对于完全支化的片段（如季碳原子或叔胺基团），如果片段包含四个非氢原子，通常会使用较小的珠子（如 S 型珠子），因为中心原子的环境暴露程度降低，其对整体相互作用的影响也相应减小 (8)。 2.3 命名约定 (Naming Conventions) Martini 3 珠子的命名遵循一套系统的规则，通常结合了其尺寸、基本化学类型、极性水平以及特殊功能：尺寸前缀: R: 代表常规尺寸 (Regular)，但在很多情况下，如果珠子名称没有明确的尺寸前缀，则默认为常规尺寸。 S: 代表小尺寸 (Small)。 T: 代表微小尺寸 (Tiny)。基本类型字母: P, N, C, Q, X, W, D，如上所述。极性/相互作用能级: 通常是一个数字（1到6，对于P, N, C, Q, X 类型），表示相对极性或相互作用强度。例如，P1 表示低极性极性珠子，P6 表示高极性极性珠子 (10)。氢键后缀: 用于描述珠子的氢键能力 (10)。 d (donor): 表示氢键给体。 a (acceptor): 表示氢键受体。 da: 表示同时具有氢键给体和受体能力。 0 (zero): 对于Q类型珠子 (如 Q0)，表示不具有特定的氢键给体或受体能力。对于P和N类型珠子，若无 ‘d’ 或 ‘a’ 后缀，通常意味着其氢键能力不是其主要特征，或作为一般极性珠子处理。电子极化效应后缀: 主要用于 C 和 X 类型珠子，以模拟邻近化学基团的诱导/共轭效应对分子片段相互作用的影响，并能捕捉优先取向和卤键能力 (8)。 e (electron-donor/enriched): 表示富电子特性。 v (electron-acceptor/vacancy): 表示缺电子特性。例如，萘中心的珠子类型为 TC5e，表示这是一个富电子的微小尺寸非极性珠子。其他特殊后缀: h: 在某些特定珠子类型中使用，例如在脂质尾链中，C5h 和 C4h 用来区分包含不同数量双键的片段 (12)。 r: 在某些溶剂模型中出现，如甲醇 (MEOH) 用 SP2r 表示 (13)。一个典型的 Martini 3 珠子名称组合了这些元素，例如：SP1d (小尺寸、极性类型、1级极性、氢键给体)，TC5e (微小尺寸、非极性类型、5级相互作用、富电子)。 2.4 Martini 3 珠子类型汇总表 (Comprehensive Table of Martini 3 Bead Types) 为了更清晰地展示 Martini 3 中常用珠子的特性，下表总结了部分代表性珠子类型及其关键属性和通常代表的化学片段。此表并非详尽无遗，更完整的列表和特定分子的参数化可以在 Martini 官方网站和相关出版物中找到 (10)。构建新分子模型时，应参考最新的官方 martini_v3.0.0.itp 文件和相关文献中的构建模块表 (8)。珠子名称 (Bead Name) 主要类型 (Main Type) 尺寸 (Size) 极性水平 (Polarity Level) 氢键 (H-bond) 其他后缀 (Other Suffix) 典型化学基团/片段 (Typical Chemical Group/Fragment) W W R N/A N/A 水 (代表4个水分子) TW W T N/A N/A 微小水 (代表2个水分子)，用于受限空间 C1 C R 1 None 饱和烷烃片段 (-CH2-CH2-CH2-CH2-) SC3 C S 3 None 脂肪环片段 (如环己烷中的 -CH2-CH2-CH2- 单元)，支链烷烃 TC5 C T 5 None 芳香环中的 -CH=CH- 片段 (如苯)，共轭体系 TC5e C T 5 None e 富电子芳香片段 (如萘的中心) P1 P R 1 Donor/Acceptor 弱极性基团，如醚 (-O-) SP2d P S 2 Donor 中等极性氢键给体，如伯醇 (-CH2OH 中的 OH 部分，若映射为S珠) TP4a P T 4 Acceptor 强极性氢键受体，如羰基 (C=O，若映射为T珠) N0 N R 0 (特殊) None 中性非极性基团，但归类于N，如某些胺的非极性部分 SN1a N S 1 Acceptor 弱中等极性氢键受体，如叔胺 (-N(CH3)-) TN4a N T 4 Acceptor 中等极性氢键受体，如醚氧 (-CH2†-O-CH2†-) Q0 Q R 0 (特殊) None 带形式电荷但无特定氢键能力的基团，或用于描述电荷离域的离子 SQd Q S (level varies) Donor 带电荷的氢键给体，如质子化的胺基 (-NH3+) TQa Q T (level varies) Acceptor 带电荷的氢键受体，如羧酸根 (-COO-) X3h X (R/S/T) 3 None h 含卤素化合物，如二氯乙烷中的 -CHCl-CH2Cl (X3h 代表一个氯原子和部分碳链) 注：上表仅为示例，具体的珠子类型选择和参数化应参考最新的 Martini 文档和相关研究。极性水平和氢键能力可能因具体的化学环境和参数化目标而有所调整。“N/A”表示不适用。 3. Martini 3 珠子参数化策略 (Martini 3 Bead Parametrization Strategy) Martini 3 珠子的参数化遵循系统性的方法，结合了“自上而下”（top-down，基于实验热力学数据）和“自下而上”（bottom-up，基于全原子模拟数据）的策略，旨在准确再现分子的物理化学性质 (2)。 3.1 非键相互作用 (Non-bonded Interactions) 非键相互作用的参数化主要目标是再现小分子在不同溶剂对之间的分配自由能 ($\Delta G_{\text{transfer}}$) (8)。常用的溶剂对包括正十六烷/水、正辛醇/水和氯仿/水等。通过拟合这些分配自由能，可以有效地校准溶质-溶剂以及溶剂-溶剂之间的交叉相互作用强度 (8)。第二个核心参数化目标是溶剂的互溶性数据，可以通过定性观察或计算混合过剩自由能来进行检验 (8)。互溶性数据同样能够反映不同种类分子间的交叉相互作用以及它们各自的自相互作用的相对强度。非键相互作用通常采用 Lennard-Jones (LJ) 势来描述： $V_{LJ}(r_{ij}) = 4 \varepsilon_{ij} \left[ \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{12} - \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{6} \right]$ 其中 rij 是珠子 i 和 j 之间的距离，$\sigma_{ij}$定义了珠子间的有效直径（相互作用为零的点），$\varepsilon_{ij}$定义了势阱深度，代表相互作用强度。对于带电荷的 Q 型和 D 型珠子，除了 LJ 相互作用外，还包含库仑相互作用： $V_C\left(r_{i j}\right)=\frac{1}{4 \pi \varepsilon_0 \varepsilon_r} \frac{q_i q_j}{r_{i j}}$ 其中$q_i$和$q_j$是珠子的电荷，$ε_0$是真空介电常数，$ε_r$是相对介电常数（在 Martini 中通常设为15，用于隐式地考虑水的屏蔽效应，但具体值可能因模拟体系而异）。 3.2 键合相互作用与几何中心映射 (Bonded Interactions and Center-of-Geometry (COG) Mapping) 与主要依赖实验数据的非键参数化不同，键合相互作用（键长、键角、二面角）的参数主要通过“自下而上”的方法获得，即参考全原子 (All-Atom, AA) 模拟得到的结构分布 (16)。 Martini 3 的一个重要改进是采用了基于“几何中心”(Center-of-Geometry, COG) 的映射规则来定义 CG 模型的键合参数，取代了 Martini 2 中常用的“质量中心”(Center of Mass, COM) 映射 (8)。COG 映射在计算分子片段中心时考虑了氢原子的位置，这使得 CG 模型能更好地保持其对应全原子参考结构的体积和形状 (8)。COM 映射有时会导致不满意的键长和过高的堆积密度，而 COG 映射则能产生更接近实际的分子性质（如溶剂可及表面积）和本体性质（如质量密度）(8)。对于接近全原子分辨率的映射（例如使用 T 型珠子），COG 映射尤为关键。例如，对于4对1映射的烷烃链，COM 和 COG 的差异不大；但对于2对1映射的苯环，两者差异显著 (8)。基于 COG 的键长可以直接从全原子模型中提取，这使得参数化过程更易于自动化。这些初始参数在需要更高精度时可以被进一步优化 (8)。常用的键合势函数包括：键长 (Bonds): 简谐势$V_b(l) = \frac{1}{2} K_b (l - l_0)^2$ 键角 (Angles): 简谐势$V_a(\theta) = \frac{1}{2} K_{\theta} (\theta - \theta_0)^2$ 二面角 (Dihedrals): 周期性势$V_d(\phi) = K_{\phi} [1 + \cos(n\phi - \phi_0)]$在某些情况下，特别是对于需要保持刚性平面结构或特定构象的分子，也会使用约束 (constraints) 或特殊势函数（如improper二面角）(8)。 4. Martini 3 映射方法学 (Martini 3 Mapping Methodology) 将全原子结构映射到粗粒化表示是构建 Martini 模型的首要步骤。Martini 3 提供了一套更一致的规则和指导原则，旨在优化 CG 模型的体积和形状表示。 4.1 基本原则 (Basic Principles) 进行原子到珠子的映射时，应遵循以下基本原则 (5)：原子分组: 通常将2-4个非氢重原子（及其相连的氢原子）映射为一个 CG 珠子。R、S、T 型珠子分别对应约4、3、2个重原子的映射。化学基团完整性: 尽量避免将特定的化学官能团（如酰胺基、羧基、完整的芳香环单元）分割到不同的珠子中，以保持其化学特性。对称性保留: 映射方案应尽可能尊重原始分子的对称性。体积与形状保持: CG 模型应能较好地再现全原子结构的体积和形状。COG 映射对此有重要贡献。珠子数量优化: 珠子的总数应被优化，目标是使每个珠子代表的重原子数与理想映射（如4:1, 3:1, 2:1）的最大偏差控制在每10个重原子中±1个非氢原子的范围内。原子共享: 在某些情况下，为了保持底层原子结构的对称性（例如在苯酚、四氢呋喃、甲苯等分子中），一个或多个原子可能被相邻的珠子共享（在一些文献的表格中用 † 标出）(8)。在从 COG 映射的全原子模拟中提取键合参数时，必须考虑到这种共享。 4.2 环状结构映射 (Mapping Ring Structures) 环状结构因其特殊的几何形状和电子特性，在 Martini 3 中有特定的映射策略，通常使用 S 型或 T 型珠子 (8)：芳香环 (Aromatic Rings): 芳香环（如苯环、萘环等）通常使用 T 型珠子进行描述，以更好地再现其平面性和堆叠行为。苯 (Benzene): 作为典型的芳香化合物，苯被模型化为三个 TC5 类型的 T 型珠子，每个珠子代表两个连续的碳原子及其相连的氢原子。TC5 是非取代芳香环中 -C=C- 基团的首选珠子类型。使用基于 COG 的键长（约 0.29 nm），可以很好地再现苯的液体密度。芳香环模型通常使用约束来连接珠子，因为其键长分布非常窄，需要非常刚性的势函数，这反过来又可能需要较短的模拟时间步长。对于更延展的刚性结构，可以考虑使用虚拟位点。脂肪环 (Aliphatic Rings): 脂肪环（如环己烷）通常使用 S 型珠子进行描述，以捕捉其相对于芳香环更大的体积感。环己烷 (Cyclohexane): 作为典型的脂肪环化合物，环己烷通常被描述为一个双 S 珠模型 (SC3-SC3)。SC3 珠子的选择基于分配数据。脂肪环模型通常使用键合相互作用（而非约束）连接，因为它们的键长分布相对较宽。使用约 0.378 nm 的键长，可以很好地再现环己烷的液体密度。 4.3 取代基与支链映射 (Mapping Substituents and Branched Chains) 对于带有取代基的环状结构或具有支链的分子，映射时需遵循以下两个主要原则 (8)：用最少数量的珠子映射所有非氢原子。尽可能保持分子的对称性、体积和形状，其中芳香环最好用 T 型珠子描述，脂肪环最好用 S 型珠子描述。例如：甲苯 (Toluene): 在苯环上增加一个甲基时，苯环原有的三个 T 型珠子中的一个会变成一个更大的 S 型珠子，以容纳额外的碳原子 (8)。乙苯 (Ethyl-benzene): 如果是乙基取代，则会为乙基额外增加一个 T 型珠子（代表乙基的两个碳原子），而苯环部分则可以保持其精确的三个 T 型珠子模型 (8)。支链烷烃/完全支化基团: 对于如新戊烷（包含5个非氢原子）这样的完全支化基团，通常会使用尺寸较小的珠子。例如，尽管有5个重原子，但由于中心碳原子被包埋，其对环境的暴露减少，因此可以使用一个 S 型珠子来代表整个新戊烷基团，或者根据具体情况进行更细致的划分 (8)。 4.4 高级模型设计策略 (Advanced Model Design Strategies) 对于具有多个芳香/脂肪环结构或复杂连接方式的小分子，Martini 3 提供了一些高级模型设计策略，常利用虚拟（相互作用）位点 (virtual sites) 来提高模型的数值稳定性和计算性能 (8)： “铰链”模型 (Hinge Model): 适用于刚性的稠合多环化合物，如萘 (Naphthalene)。一个简单的由5个 TC5 珠子通过约束连接的萘模型在凝聚相中可能导致数值不稳定。 “铰链”结构使用4个外部珠子，并将中心的一个或多个珠子描述为虚拟相互作用位点（其位置由构建粒子定义，受力会分配给构建粒子，质量均匀分配给构建粒子）。这种方法减少了约束数量，提高了数值稳定性和模拟速度。通常还会施加一个不当二面角来保持铰链模型的平面性。 “分而治之”模型 (Divide and Conquer): 适用于由刚性平面片段组成的任意长链，且需要控制片段间的相对二面角，这在小分子和共轭聚合物（如2,2’-联噻吩）中很常见。例如，两个噻吩环各用三个 T 型珠子描述。为了连接它们并控制二面角，可以在每个噻吩环的几何中心使用两个虚拟的非相互作用的哑位点 (dummy sites)，并通过简谐键连接这两个哑位点。然后可以在这些哑位点和每个噻吩环上的两个粒子（如硫原子）之间施加二面角势。 “分子转角”模型 (Molecular Turn): 用于处理通过 sp2 杂化碳连接的环系统，这种连接方式会产生一个“分子转角”（如某些具有特定扭转行为的分子）。为了保持扭转运动的正确旋转轴，需要特别注意。通常会使用虚拟哑位点：每个环体系的 COG 处各一个，第三个位于连接的 sp2 杂化碳上。通过在这些虚拟位点之间施加简谐键和角势，并辅以适当的二面角势和不当二面角势来维持正确的几何构型和旋转自由度。这些高级策略体现了 Martini 3 在处理复杂分子结构方面的灵活性和精确性。 5. 总结与拓展资源 (Conclusion and Further Resources) 5.1 总结 (Summary) Martini 3 力场通过引入新的珠子类型、更精细的尺寸划分 (R, S, T)、系统的命名规则（包含极性、氢键能力、电子特性等后缀）以及改进的参数化策略（特别是基于几何中心 COG 的映射），显著提升了粗粒化模拟的准确性和适用范围 (1)。其核心优势在于能够在保持较高计算效率的同时，捕捉到关键的化学物理特性，从而能够模拟更大尺度和更长时间尺度的分子过程。在对新分子（尤其是如本教程中所示的复杂聚合物）进行 Martini 3 映射时，关键步骤包括：仔细分析全原子化学结构，识别关键官能团。遵循原子分组（2-4个重原子/珠子）、化学基团完整性、对称性和体积/形状保持等基本映射原则。参考 Martini 3 珠子类型表和命名规则，为每个分子片段选择最合适的珠子类型和尺寸。对于缺乏直接预参数化珠子的特殊基团（如本例中的N-氧化物），需基于其化学物理特性（极性、氢键、分配行为等）类比选择最接近的现有珠子，或进行审慎的重新参数化。定义珠子间的键合连接。通过这些步骤，可以为目标分子构建合理的 Martini 3 粗粒化模型，为后续的分子动力学模拟打下坚实基础。 5.2 拓展资源 (Further Resources) 为了更深入地学习和应用 Martini 3 力场，以下资源非常宝贵： Martini 官方网站: http://cgmartini.nl (2)。这里可以找到最新的力场文件、教程、FAQ 以及已参数化的分子拓扑数据库（包括脂质、蛋白质、糖类、溶剂和小分子等 (8)）。主要出版物: Souza, P.C.T., Alessandri, R., Barnoud, J. et al. Martini 3: a general purpose force field for coarse-grained molecular dynamics. Nat Methods 18, 382–388 (2021). (1) (Martini 3 的奠基性论文)。 Alessandri, R., Souza, P.C.T., Thallmair, S. et al. A coarse-grained force field for small molecules: Martini 3. ChemRxiv (2021). (8) (针对小分子参数化的重要参考，包含大量构建模块信息)。模拟软件包: Martini 力场广泛应用于 GROMACS (4)。NAMD 等其他软件包也有相应的支持或转换工具 (26)。辅助工具: Martinize (或 Martinize2): 用于将全原子蛋白质（以及其他分子）结构转换为 Martini CG 模型的常用脚本 (5)。 Insane.py: 用于快速搭建复杂膜体系的脚本 (24)。 Polyply: 用于生成聚合物拓扑的工具 (29)。 MartiniGlass: 用于 VMD 中可视化 Martini 模型的 Python 包 (23)。力场参数下载: Martini 3 核心参数文件 (martini_v300.zip 或类似名称) 可从官方网站下载，其中包含了珠子定义 (.itp 文件)、相互作用矩阵以及多种已参数化分子的拓扑文件 (14)。小分子数据库通常托管在 GitHub 等代码仓库中，如 ricalessandri/Martini3-small-molecules (8)。 Marrink实验室的 GitHub 仓库 (marrink-lab/martini-forcefields) 也是获取最新参数和分子拓扑的重要来源 (29)。利用这些资源，研究者可以有效地将 Martini 3 应用于广泛的化学和生物物理问题研究中，探索复杂体系的结构、动态和热力学性质。

Molecular Dynamics · 2025-11-02

Mendelevium

Contact

FF & Algorithm